\lstset{language=Octave, commentstyle=\color{gray}, stringstyle=\color{red!60!black}}
\chapter{FDTD Code for 1D Wave Propagation in Free Space}
\label{App:1DFDTDFreeSpaceMatlab}
\begin{lstlisting}
% Simulation parameters.
SIZE = 512; % No. of spatial steps
MaxTime = 512; % No. of time steps
PulseWidth = 128; % Width of Gaussian Pulse
imp0 = 377.0; % Impedence of free space

% Constants.
c = 3e8;
pi = 3.141592654;
e0 = (1e-9)/(36*pi);
u0 = (1e-7)*4*pi;

dt = 1e-12;
dz = 3e-4;
Sc = c * dt/dz

% Initialization.
Ex = zeros ( SIZE, MaxTime ); % z-component of E-field
Hy = zeros ( SIZE, MaxTime ); % y-component of H-field
PLOT1(1) = 0; % Data for plotting.

% Outer loop for time-stepping.
tic
for q = 2:MaxTime
    % Calculation of Hy using update difference equation for Hy. This is time step q.
    Hy(1:SIZE-1,q) = Hy(1:SIZE-1,q-1) + ( ( Ex(1:SIZE-1,q-1) - Ex(2:SIZE,q-1) ) * dt/(u0*dz) );
    % ABC for H at SIZE.
    Hy(SIZE,q) = Hy(SIZE-1,q-1) + (Sc-1)/(Sc+1)*(Hy(SIZE-1,q) - Hy(SIZE,q-1) );
    
    % Calculation of Ez using updated difference equation for Ez. This is time step q+1/2.
    Ex(2:SIZE,q) = Ex(2:SIZE, q-1) + ( dt/(e0*dz)*(Hy(1:SIZE-1, q) - Hy(2:SIZE, q)) );
    % ABC for E at 1.
    Ex(1,q) = Ex(2,q-1) + (Sc-1)/(Sc+1)*(Ex(2,q) - Ex(2,q-1));
    
    % Activating a plane-wave source.
    Ex(100,q) = Ex(100,q) + exp( -1*(q-PulseWidth)^2/PulseWidth ) * Sc;
end
toc

% Simulation animation.
for i=1:MaxTime
    figure (2)
    plot ( Ex(:,i) )
end

\end{lstlisting}

\chapter{1D FDTD Demonstrating PBC}
\label{App:1DFDTDFreeSpacePBCMatlab}
\begin{lstlisting}
% Simulation parameters.
SIZE = 1024; % No. of spatial steps
MaxTime = 2024; % No. of time steps
imp0 = 377.0; % Impedence of free space

% Constants.
c = 3e8;
pi = 3.141592654;
e0 = (1e-9)/(36*pi);
u0 = (1e-7)*4*pi;

dt = 0.8e-12;
dx = 3e-4;
Sc = c * dt/dx

% Initialization.
Ez = zeros ( SIZE, MaxTime ); % z-component of E-field
Hy = zeros ( SIZE, MaxTime ); % y-component of H-field
PLOT1(1) = 0; % Data for plotting.

% Outer loop for time-stepping.
tic
for q = 2:MaxTime
    % Calculation of Hy using update difference equation for Hy. This is time step q.
    Hy(1:SIZE-1,q) = Hy(1:SIZE-1,q-1) + ( ( Ez(2:SIZE,q-1) - Ez(1:SIZE-1,q-1) ) * dt/(u0*dx) );
    % PBC for H at SIZE.
    Hy(SIZE,q) = Hy(SIZE,q-1) + ( ( Ez(1,q-1) - Ez(SIZE,q-1) ) * dt/(u0*dx) );
    
    % Calculation of Ez using updated difference equation for Ez. This is time step q+1/2.
    Ez(2:SIZE,q) = Ez(2:SIZE, q-1) + ( dt/(e0*dx)*(Hy(2:SIZE, q) - Hy(1:SIZE-1, q)) );
    % PBC for E at 1.
    Ez(1,q) = Ez(1,q-1) + ( dt/(e0*dx)*(Hy(1, q) - Hy(SIZE, q)) );
    
    % Activating a plane-wave source.
    Ez(100,q) = Ez(100,q) + exp( -1*((q-256)^2)/1000 ) * Sc;
end
toc

% Simulation animation.
for i=1:MaxTime
    figure (2)
    plot ( Ez(:,i) )
end

\end{lstlisting}
\chapter{1D DNG Simulation}
\label{App:1DFDTDDNGMatlab}
\begin{lstlisting}
clc
clear all

% Simulation parameters.
SIZE = 4*1024; % No. of spatial steps
SlabLeft = round(SIZE/3); % Location of left end of Slab.
SlabRight = round(2*SIZE/3); % Location of right end of Slab
MaxTime = 4*SIZE; % No. of time steps
PulseWidth = round(SIZE/8); % Controls width of Gaussian Pulse
td = PulseWidth; % Temporal delay in pulse.
source = 10; % Location of source
SnapshotInterval = 16; % Amount of time delay between snaps.

% Constants.
c = 3e8;
pi = 3.141592654;
e0 = (1e-9)/(36*pi);
u0 = (1e-7)*4*pi;
dt = 0.5e-11;
dz = 3e-3;
Sc = c * dt/dz
l = PulseWidth*dz;
f = c/l
fmax = 1/(2*dt)
w = 2*pi*f;
k0 = w/c; % Free space wave number.

% Choice of source.
% 1. Gaussian 2. Sine wave 3. Ricker wavelet
SourceChoice = 1;
if SourceChoice == 3
    fp = f; % Peak frequency
    dr = PulseWidth*dt*2; % Delay
end

% Initialization.
Ex = zeros(SIZE, 3); % x-component of E-field
Dx = zeros(SIZE, 3); % x-component of D
Hy = zeros(SIZE, 3); % y-component of H-field
By = zeros(SIZE, 3); % y-component of B

% Incident and Transmitted Fields.
Exi = zeros(MaxTime, 1);
Ext = zeros(MaxTime, 1);
Extt = zeros(MaxTime, 1);
x1 = SlabLeft+1; % Position of observation.

% Refractive Index calculations.
Z1 = SlabLeft + 50;
z1 = Z1*dz;
Z2 = SlabLeft + 60;
z2 = Z2*dz;
Exz1 = zeros(MaxTime, 1);
Exz2 = zeros(MaxTime, 1);

einf = ones(SIZE,1);
einf(SlabLeft:SlabRight) = 1; % einf(Drude) or er in slab.
uinf = ones(SIZE,1);
uinf(SlabLeft:SlabRight) = 1; % uinf(Drude) or ur in slab.
wpesq = zeros(SIZE,1);
wpesq(SlabLeft:SlabRight) = 2*w^2; % DNG(Drude) value of wpe squared in slab.
wpmsq = zeros(SIZE,1);
wpmsq(SlabLeft:SlabRight) = 2*w^2; % DNG(Drude) value of wpm squared in slab.
ge = zeros(SIZE,1);
ge(SlabLeft:SlabRight) = w/16; % Electric collision frequency in slab.
gm = zeros(SIZE,1);
gm(SlabLeft:SlabRight) = w/16; % Magnetic collision frequency in slab.

a0 = (4*dt^2)./(e0*(4*einf+dt^2*wpesq+2*dt*einf.*ge));
a = (1/dt^2)*a0;
b = (1/(2*dt))*ge.*a0;
c = (e0/dt^2)*einf.*a0;
d = (-1*e0/4).*wpesq.*a0;
e = (1/(2*dt))*e0*einf.*ge.*a0;
am0 = (4*dt^2)./(u0*(4*uinf+dt^2*wpmsq+2*dt*uinf.*gm));
am = (1/dt^2)*am0;
bm = (1/(2*dt))*gm.*am0;
cm = (u0/dt^2)*uinf.*am0;
dm = (-1*u0/4).*wpmsq.*am0;
em = (1/(2*dt))*u0*uinf.*gm.*am0;

ExSnapshots = zeros(SIZE, MaxTime/SnapshotInterval); % Data for plotting.
frame = 1;
n1 = 1;
n2 = 2;
% Outer loop for time-stepping.
tic
% Test loop for incident field in free space.
for q = 0:MaxTime
   
    % Calculation of Hy using update difference equation for Hy. This is time step q.
    Hy(1:SIZE-1,n2) = Hy(1:SIZE-1,n1) + ( ( Ex(1:SIZE-1,n1) - Ex(2:SIZE,n1) ) * dt/(u0*dz) );
        
    % ABC for H at SIZE.
    Hy(SIZE,n2) = Hy(SIZE-1,n1) + (Sc-1)/(Sc+1)*(Hy(SIZE-1,n2) - Hy(SIZE,n1) );

    % Calculation of Ex using updated difference equation for Ex. This is time step q+1/2.
    Ex(2:SIZE,n2) = Ex(2:SIZE, n1) + ( dt/(e0*dz)*(Hy(1:SIZE-1, n2) - Hy(2:SIZE, n2)) );
    
    % ABC for E at 1.
    Ex(1,n2) = Ex(2,n1) + (Sc-1)/(Sc+1)*(Ex(2,n2) - Ex(2,n1));
    
    % Source.
    if SourceChoice == 1
    Ex(source,n2) = Ex(source,n2) + exp( -1*((q-td)/(PulseWidth/4))^2 ) * Sc;
    elseif SourceChoice == 2
    Ex(source,n2) = Ex(source,n2) + sin(2*pi*f*(q)*dt) * Sc;
    elseif SourceChoice == 3
    Ex(source,n2) = Ex(source,n2) + (1-2*(pi*fp*(q*dt-dr))^2)*exp(-1*(pi*fp*(q*dt-dr))^2) * Sc;
    end
    
    Exi(q+1) = Ex(x1,n2); % Incident field is left of slab.
    
    temp = n1;
    n1 = n2;
    n2 = temp;
end

% Reinitialization of fields for actual simulation.
Ex = zeros(SIZE, 3); % x-component of E-field
Hy = zeros(SIZE, 3); % y-component of H-field
% Actual simulation with scatterer.
for q = 0:MaxTime
    
    % Storing past fields.
    Ex(:,3) = Ex(:,n2);
    Dx(:,3) = Dx(:,n2);
    Hy(:,3) = Hy(:,n2);
    By(:,3) = By(:,n2);
    
    % Calculation of Hy using update difference equation for Hy. This is time step q.
    By(1:SIZE-1,n2) = By(1:SIZE-1,n1) + ( ( Ex(1:SIZE-1,n1) - Ex(2:SIZE,n1) ) * dt/(dz) );
    Hy(:,n2) = am.*(By(:,n2)-2*By(:,n1)+By(:,3))+bm.*(By(:,n2)-By(:,3))+cm.*(2*Hy(:,n1)-Hy(:,3))+dm.*(2*Hy(:,n1)+Hy(:,3))+em.*(Hy(:,3));
    
    % ABC for H at SIZE.
    Hy(SIZE,n2) = Hy(SIZE-1,n1) + (Sc-1)/(Sc+1)*(Hy(SIZE-1,n2) - Hy(SIZE,n1) );
    By(SIZE,n2) = u0*Hy(SIZE,n2);

    % Calculation of Ex using updated difference equation for Ex. This is time step q+1/2.
    Dx(2:SIZE,n2) = Dx(2:SIZE, n1) + ( dt/(dz)*(Hy(1:SIZE-1, n2) - Hy(2:SIZE, n2)) );
    Ex(:,n2) = a.*(Dx(:,n2)-2*Dx(:,n1)+Dx(:,3))+b.*(Dx(:,n2)-Dx(:,3))+c.*(2*Ex(:,n1)-Ex(:,3))+d.*(2*Ex(:,n1)+Ex(:,3))+e.*(Ex(:,3));
    
    % ABC for E at 1.
    Ex(1,n2) = Ex(2,n1) + (Sc-1)/(Sc+1)*(Ex(2,n2) - Ex(2,n1));
    Dx(1,n2) = e0*Ex(1,n2);
    
    % Source.
    if SourceChoice == 1
        Ex(source,n2) = Ex(source,n2) + exp( -1*((q-td)/(PulseWidth/4))^2 ) * Sc;
    elseif SourceChoice == 2
        Ex(source,n2) = Ex(source,n2) + sin(2*pi*f*(q)*dt) * Sc;
    elseif SourceChoice == 3
        Ex(source,n2) = Ex(source,n2) + (1-2*(pi*fp*(q*dt-dr))^2)*exp(-1*(pi*fp*(q*dt-dr))^2) * Sc;
    end
    Dx(source,n2) = e0*Ex(source,n2);

    if mod(q,SnapshotInterval) == 0
        ExSnapshots(:,frame) = Ex(:,n2);
        frame=frame+1;
    end
    
    Ext(q+1) = Ex(x1,n2);
    Extt(q+1) = Ex(SlabRight+10,n2);
    
    % Fields for calculation of refractive index.
    Exz1(q+1) = Ex(Z1, n2);
    Exz2(q+1) = Ex(Z2, n2);
    
    temp = n1;
    n1 = n2;
    n2 = temp;
end
toc

% Postprocessing.
Fs = 1/dt;                    % Sampling frequency
T = dt;                       % Sample time
L = length(Exi);              % Length of signal
t = (0:L-1)*T;                % Time vector
fspan = 100;                  % Points to plot in frequency domain

figure(1)
subplot(211)
plot(Fs*t, Exi, 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Incident Field', 'FontSize', 12, 'FontWeight', 'b')
xlabel('time', 'FontSize', 11, 'FontWeight', 'b')
grid on
figure(2)
subplot(211)
plot(Fs*t, Ext, 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Transmitted Field', 'FontSize', 12, 'FontWeight', 'b')
xlabel('time', 'FontSize', 11, 'FontWeight', 'b')
grid on
figure(3)
subplot(211)
plot(Fs*t, Extt, 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Transmitted Field Beyond Slab', 'FontSize', 12, 'FontWeight', 'b')
xlabel('time', 'FontSize', 11, 'FontWeight', 'b')
grid on

NFFT = 2^nextpow2(L); % Next power of 2 from length of Exi
% Incident and Transmitted fields.
EXI = fft(Exi,NFFT)/L;
EXT = fft(Ext,NFFT)/L;
EXTT = fft(Extt,NFFT)/L;
% Refractive index calculations.
EXZ1 = fft(Exz1,NFFT)/L;
EXZ2 = fft(Exz2,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure(1)
subplot(212)
EXIp = 2*abs(EXI(1:NFFT/2+1));
plot(f(1:fspan), EXIp(1:fspan), 'LineWidth', 2.0, 'Color', 'r')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Single-Sided Amplitude Spectrum of Exi(t)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frequency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('|EXI(f)|', 'FontSize', 11, 'FontWeight', 'b')
grid on
figure(2)
subplot(212)
EXTp = 2*abs(EXT(1:NFFT/2+1));
plot(f(1:fspan), EXTp(1:fspan), 'LineWidth', 2.0, 'Color', 'r')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Single-Sided Amplitude Spectrum of Ext(t)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frequency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('|EXT(f)|', 'FontSize', 11, 'FontWeight', 'b')
grid on

figure(3)
subplot(212)
EXTTp = 2*abs(EXTT(1:NFFT/2+1));
plot(f(1:fspan), EXTTp(1:fspan), 'LineWidth', 2.0, 'Color', 'r')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Single-Sided Amplitude Spectrum of Extt(t)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frequency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('|EXT(f)|', 'FontSize', 11, 'FontWeight', 'b')
grid on

% Transmission Coefficient.
figure(4)
subplot(211)
TAU = abs(EXT(1:NFFT/2+1)./EXI(1:NFFT/2+1));
plot(f(1:fspan), TAU(1:fspan), 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Transmission Coefficient', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frequency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('|EXT(f)/EXI(f)|', 'FontSize', 11, 'FontWeight', 'b')
axis([-1 1 -2 2])
axis 'auto x'
grid on
subplot(212)
plot(f(1:fspan), 1-TAU(1:fspan), 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Reflection Coefficient', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frequency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('1-|EXT(f)/EXI(f)|', 'FontSize', 11, 'FontWeight', 'b')
axis([-1 1 -2 2])
axis 'auto x'
grid on

% Refractive Index calculations.
nFDTD = (1/(1i*k0*(z1-z2))).*log(EXZ2(1:NFFT/2+1)./EXZ1(1:NFFT/2+1));
figure(5)
subplot(211)
plot(f(1:fspan), real(nFDTD(1:fspan)), 'LineWidth', 2.0, 'Color', 'b');
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Refractive index re(n)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frequency (Hz)', 'FontSize', 11)
ylabel('re(n)', 'FontSize', 11)
grid on
subplot(212)
plot(f(1:fspan), imag(nFDTD(1:fspan)), 'LineWidth', 2.0, 'Color', 'r');
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Refractive index im(n)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frequency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('im(n)', 'FontSize', 11, 'FontWeight', 'b')
grid on

% Simulation animation.
for i=1:frame-1
    figure (6)
    % Scatterer boundaries.
    hold off
    plot([SlabLeft SlabLeft], [-1 1], 'Color', 'r');
    hold on
    plot([SlabRight SlabRight], [-1 1], 'Color', 'r');
    plot(ExSnapshots(:,i), 'LineWidth', 2.0, 'Color', 'b');
    set(gca, 'FontSize', 10, 'FontWeight', 'b')
    axis([0 SIZE -1 1])
    title('Time Domain Simulation', 'FontSize', 12, 'FontWeight', 'b')
    xlabel('Spatial step (k)', 'FontSize', 11, 'FontWeight', 'b')
    ylabel('Electric field (Ex)', 'FontSize', 11, 'FontWeight', 'b')
    grid on
end
\end{lstlisting}
\chapter{2D DNG Simulation}
\label{App:2DFDTDDNGMatlab}
\begin{lstlisting}
clc
clear all

% Simulation parameters.
SizeI = 512; % No. of spatial steps in x direction.
SizeJ = 512; % No. of spatial steps in y direction.
PMLw = 50; % Width of PML layer.
SlabLeft = round(SizeJ/3+PMLw); % Location of left end of Slab.
SlabRight = round(2*SizeJ/3+PMLw); % Location of right end of Slab
MaxTime = 4*SizeJ; % No. of time steps
PulseWidth = round(SizeJ/8); % Controls width of Gaussian Pulse
td = PulseWidth; % Temporal delay in pulse.
SnapshotResolution = 1; % Snapshot resolution. 1 is best.
SnapshotInterval = 4; % Amount of time delay between snaps.
% Choice of source.
% 1. Gaussian 2. Sine wave 3. Ricker wavelet
SourceChoice = 1;
SourcePlane = 1; % Is the source a plane wave. 0. = Omni 1. Plane-wave.
SourceLocationX = SizeI/2; % X Location of source. Only used for an omni-source.
SourceLocationY = PMLw+10; % Y Location of source.

% Constants.
c = 3e8;
pi = 3.141592654;
e0 = (1e-9)/(36*pi);
u0 = (1e-7)*4*pi;

dt = 0.5e-11;
delta = 3e-3;
Sc = c * dt/delta

l = PulseWidth*delta;
f = c/l
fmax = 1/(2*dt)
w = 2*pi*f;
k0 = w/c; % Free space wave number.
% Ricker wavelet parameters.
if SourceChoice == 3
    fp = f; % Peak frenuency
    dr = PulseWidth*dt*2; % Delay
end
% PML arrays.
PsiEzX = zeros(SizeI, SizeJ+2*PMLw);
PsiEzY = zeros(SizeI, SizeJ+2*PMLw);
PsiHyX = zeros(SizeI, SizeJ+2*PMLw);
PsiHxY = zeros(SizeI, SizeJ+2*PMLw+1);

% PML parameters.
kapp = 1;
a = 0.0004;
sig = 0.045;
% Electric.
kappex = kapp;
kappey = kapp;
aex = a;
aey = a;
sigex = 0;
sigey = sig;
bex = exp(-1*(aex/e0+sigex/(kappex*e0))*dt);
bey = exp(-1*(aey/e0+sigey/(kappey*e0))*dt);
Cex = (bex-1)*sigex/(sigex*kappex+kappex^2*aex);
Cey = (bey-1)*sigey/(sigey*kappey+kappey^2*aey);
% Magnetic.
kappmx = kapp;
kappmy = kapp;
amx = a;
amy = a;
sigmx = 0;
sigmy = u0/e0*sig;
bmx = exp(-1*(amx/u0+sigmx/(kappmx*u0))*dt);
bmy = exp(-1*(amy/u0+sigmy/(kappmy*u0))*dt);
Cmx = (bmx-1)*sigmx/(sigmx*kappmx+kappmx^2*amx);
Cmy = (bmy-1)*sigmy/(sigmy*kappmy+kappmy^2*amy);

% Initialization.
Ez = zeros(SizeI, SizeJ+2*PMLw, 3); % z-component of E-field
Dz = zeros(SizeI, SizeJ+2*PMLw, 3); % z-component of D
Hx = zeros(SizeI, SizeJ+2*PMLw+1, 3); % x-component of H-field
Bx = zeros(SizeI, SizeJ+2*PMLw+1, 3); % x-component of B
Hy = zeros(SizeI, SizeJ+2*PMLw, 3); % y-component of H-field
By = zeros(SizeI, SizeJ+2*PMLw, 3); % y-component of B

% Incident and Transmitted Fields.
Ezi = zeros(MaxTime, 1);
Ezt = zeros(MaxTime, 1);
Eztt = zeros(MaxTime, 1);
x1 = SlabLeft+1; % Position of observation.

% Refractive Index calculations.
Y1 = SlabLeft + 5;
y1 = Y1*delta;
Y2 = SlabLeft + 6;
y2 = Y2*delta;
Ezy1 = zeros(MaxTime, 1);
Ezy2 = zeros(MaxTime, 1);

einf = ones(SizeI,SizeJ+2*PMLw+1);
einf(:,SlabLeft:SlabRight) = 1; % einf(Drude) or er in slab.
uinf = ones(SizeI,SizeJ+2*PMLw+1);
uinf(:,SlabLeft:SlabRight) = 1; % uinf(Drude) or ur in slab.
wpesn = zeros(SizeI,SizeJ+2*PMLw+1);
wpesn(:,SlabLeft:SlabRight) = 2*w^2; % DNG(Drude) value of wpe snuared in slab.
wpmsn = zeros(SizeI,SizeJ+2*PMLw+1);
wpmsn(:,SlabLeft:SlabRight) = 2*w^2; % DNG(Drude) value of wpm snuared in slab.
ge = zeros(SizeI,SizeJ+2*PMLw+1);
ge(:,SlabLeft:SlabRight) = w/32; % Electric collision frenuency in slab.
gm = zeros(SizeI,SizeJ+2*PMLw+1);
gm(:,SlabLeft:SlabRight) = w/32; % Magnetic collision frenuency in slab.

ae0 = (4*dt^2)./(e0*(4*einf+dt^2*wpesn+2*dt*einf.*ge));
ae = (1/dt^2)*ae0;
be = (1/(2*dt))*ge.*ae0;
ce = (e0/dt^2)*einf.*ae0;
de = (-1*e0/4).*wpesn.*ae0;
ee = (1/(2*dt))*e0*einf.*ge.*ae0;
am0 = (4*dt^2)./(u0*(4*uinf+dt^2*wpmsn+2*dt*uinf.*gm));
am = (1/dt^2)*am0;
bm = (1/(2*dt))*gm.*am0;
cm = (u0/dt^2)*uinf.*am0;
dm = (-1*u0/4).*wpmsn.*am0;
em = (1/(2*dt))*u0*uinf.*gm.*am0;

EzSnapshots = zeros(SizeI/SnapshotResolution, (SizeJ+2*PMLw)/SnapshotResolution, MaxTime/SnapshotInterval); % Data for plotting.
frame = 1;

np = 1;
n0 = 2;
nf = 3;
linecount = 0;
% Outer loop for time-stepping.
fprintf ( 1, '\rDry run started! \n');
tic
% Test loop for incident field in free space.
for n = 0:MaxTime

    % Progress indicator.
    if mod(n,2) == 0
        fprintf(1, repmat('\b',1,linecount));
        linecount = fprintf(1, '%g %%', (n*100)/MaxTime );
    end

    % ========================= Bx and Hx =============================
    % Hx Psi array.
    x=1:SizeI;
    y=2:SizeJ+2*PMLw;
    PsiHxY(x,y) = (Cmy/delta)*(-Ez(x,y,n0) + Ez(x,y-1,n0)) + bmy*PsiHxY(x,y);
    % Bx in normal space.
    y=(2+PMLw):((SizeJ+2*PMLw+1)-PMLw-1);
    Bx(x,y,nf) = Bx(x,y,n0) + (-Ez(x,y,n0) + Ez(x,y-1,n0)) * dt/delta;
    if PMLw > 0
        % Bx in lower PML layer.
        y=2:PMLw+1;
        Bx(x,y,nf) = Bx(x,y,n0) + dt*((1/kappmy)*(-Ez(x,y,n0) + Ez(x,y-1,n0)) * 1/delta + PsiHxY(x,y));
        % Bx in upper PML layer.
        y=(SizeJ+2*PMLw+1)-PMLw:(SizeJ+2*PMLw);
        Bx(x,y,nf) = Bx(x,y,n0) + dt*((1/kappmy)*(-Ez(x,y,n0) + Ez(x,y-1,n0)) * 1/delta + PsiHxY(x,y));
    end
    Hx(:,:,nf) = Bx(:,:,nf)./(u0*uinf);

    % ========================= By and Hy =============================
    % Hy Psi array.
    x=1:SizeI-1;
    y=1:SizeJ+2*PMLw;
    PsiHyX(x,y) = (Cmx/delta)*(Ez(x+1,y,n0)-Ez(x,y,n0)) + bmx*PsiHyX(x,y);
    PsiHyX(SizeI,y) = (Cmx/delta)*(Ez(1,y,n0)-Ez(SizeI,y,n0)) + bmx*PsiHyX(SizeI,y);
    % By in normal space.
    y=(1+PMLw):(SizeJ+2*PMLw)-PMLw;
    By(x,y,nf) = By(x,y,n0) + (Ez(x+1,y,n0) - Ez(x,y,n0)) * dt/delta;
    By(SizeI,y,nf) = By(SizeI,y,n0) + (Ez(1,y,n0) - Ez(SizeI,y,n0)) * dt/delta; % PBC
    if PMLw > 0
        % By in lower PML layer.
        y=1:PMLw;
        By(x,y,nf) = By(x,y,n0) + dt*((1/kappmx)*(Ez(x+1,y,n0) - Ez(x,y,n0)) * 1/delta + PsiHyX(x,y));
        By(SizeI,y,nf) = By(SizeI,y,n0) + dt*((1/kappmx)*(Ez(1,y,n0) - Ez(SizeI,y,n0)) * 1/delta + PsiHyX(SizeI,y)); % PBC
        % By in upper PML layer.
        y=(SizeJ+2*PMLw)-PMLw+1:(SizeJ+2*PMLw);
        By(x,y,nf) = By(x,y,n0) + dt*((1/kappmx)*(Ez(x+1,y,n0) - Ez(x,y,n0)) * 1/delta + PsiHyX(x,y));
        By(SizeI,y,nf) = By(SizeI,y,n0) + dt*((1/kappmx)*(Ez(1,y,n0) - Ez(SizeI,y,n0)) * 1/delta + PsiHyX(SizeI,y)); % PBC
    end
    Hy(:,:,nf) = By(:,:,nf)./(u0*uinf(:,1:SizeJ+2*PMLw));

    % ========================= Dz and Ez =============================
    % Psi arrays.
    x=2:SizeI;
    y=1:SizeJ+2*PMLw;
    PsiEzX(x,y) = (Cex/delta)*(Hy(x,y,nf)-Hy(x-1,y,nf)) + bex*PsiEzX(x,y);
    PsiEzX(1,y) = (Cex/delta)*(Hy(1,y,nf)-Hy(SizeI,y,nf)) + bex*PsiEzX(1,y); % PBC
    PsiEzY(x,y) = (Cey/delta)*(-Hx(x,y+1,nf)+Hx(x,y,nf)) + bey*PsiEzY(x,y);
    PsiEzY(1,y) = (Cey/delta)*(-Hx(1,y+1,nf)+Hx(1,y,nf)) + bey*PsiEzY(1,y); % PBC
    % Dz in Normal Space.
    y=(1+PMLw):((SizeJ+2*PMLw)-PMLw);
    Dz(x,y,nf) = Dz(x,y,n0) + (Hy(x,y,nf)-Hy(x-1,y,nf)-Hx(x,y+1,nf)+Hx(x,y,nf)) * dt/delta;
    Dz(1,y,nf) = Dz(1,y,n0) + (Hy(1,y,nf)-Hy(SizeI,y,nf)-Hx(1,y+1,nf)+Hx(1,y,nf)) * dt/delta; % PBC
    if PMLw > 0
        % Dz in lower PML layer.
        y=1:PMLw;
        Dz(x,y,nf) = Dz(x,y,n0) + dt*(((1/kappex)*(Hy(x,y,nf)-Hy(x-1,y,nf))+(1/kappey)*(-Hx(x,y+1,nf)+Hx(x,y,nf))) * 1/delta + PsiEzX(x,y) + PsiEzY(x,y));
        Dz(1,y,nf) = Dz(1,y,n0) + dt*(((1/kappex)*(Hy(1,y,nf)-Hy(SizeI,y,nf))+(1/kappey)*(-Hx(1,y+1,nf)+Hx(1,y,nf))) * 1/delta + PsiEzX(1,y) + PsiEzY(1,y)); % PBC
        % Dz in upper PML layer.
        y=(SizeJ+2*PMLw)-PMLw+1:(SizeJ+2*PMLw);
        Dz(x,y,nf) = Dz(x,y,n0) + dt*(((1/kappex)*(Hy(x,y,nf)-Hy(x-1,y,nf))+(1/kappey)*(-Hx(x,y+1,nf)+Hx(x,y,nf))) * 1/delta + PsiEzX(x,y) + PsiEzY(x,y));
        Dz(1,y,nf) = Dz(1,y,n0) + dt*(((1/kappex)*(Hy(1,y,nf)-Hy(SizeI,y,nf))+(1/kappey)*(-Hx(1,y+1,nf)+Hx(1,y,nf))) * 1/delta + PsiEzX(1,y) + PsiEzY(1,y)); % PBC
    end
    Ez(:,:,nf) = Dz(:,:,nf)./(e0*einf(:,1:SizeJ+2*PMLw));

    % ====================== Source ===================
    if SourcePlane == 1
        x = 1:SizeI;
        y = SourceLocationY;        
    else
        x = SourceLocationX;
        y = SourceLocationY;
    end
    % Source.
    if SourceChoice == 1
        Ez(x,y,nf) = Ez(x,y,nf) + exp( -1*((n-td)/(PulseWidth/4))^2 ) * Sc;
    elseif SourceChoice == 2
        Ez(x,y,nf) = Ez(x,y,nf) + sin(2*pi*f*(n)*dt) * Sc;
    elseif SourceChoice == 3
        Ez(x,y,nf) = Ez(x,y,nf) + (1-2*(pi*fp*(n*dt-dr))^2)*exp(-1*(pi*fp*(n*dt-dr))^2) * Sc;
    end
    Dz(x,y,nf) = e0*Ez(x,y,nf);

    Ezi(n+1) = Ez(SizeI/2,SlabLeft+1,nf); % Incident field is left of slab.
    
    if (mod(n, SnapshotInterval) == 0)
        EzSnapshots(:,:,n/SnapshotInterval+1) = Ez(1+(0:SnapshotResolution:(SizeI-1)), 1+(0:SnapshotResolution:((SizeJ+2*PMLw)-1)), nf);
    end

    np = mod(np, 3)+1;
    n0 = mod(n0, 3)+1;
    nf = mod(nf, 3)+1;
end
fprintf ( 1, '\rDry run complete! \n');
toc
% Reinitialization of fields for actual simulation.
Ez = zeros(SizeI, SizeJ+2*PMLw, 3); % z-component of E-field
Dz = zeros(SizeI, SizeJ+2*PMLw, 3); % z-component of D
Hx = zeros(SizeI, SizeJ+2*PMLw+1, 3); % x-component of H-field
Bx = zeros(SizeI, SizeJ+2*PMLw+1, 3); % x-component of B
Hy = zeros(SizeI, SizeJ+2*PMLw, 3); % y-component of H-field
By = zeros(SizeI, SizeJ+2*PMLw, 3); % y-component of B

% Actual simulation with scatterer.
fprintf ( 1, 'Simulation started... \n');
np = 1;
n0 = 2;
nf = 3;
linecount = 0;
% Outer loop for time-stepping.
tic
% Test loop for incident field in free space.
for n = 0:MaxTime

    % Progress indicator.
    if mod(n,2) == 0
        fprintf(1, repmat('\b',1,linecount));
        linecount = fprintf(1, '%g %%', (n*100)/MaxTime );
    end

    % ========================= Bx and Hx =============================
    % Hx Psi array.
    x=1:SizeI;
    y=2:SizeJ+2*PMLw;
    PsiHxY(x,y) = (Cmy/delta)*(-Ez(x,y,n0) + Ez(x,y-1,n0)) + bmy*PsiHxY(x,y);
    % Bx in normal space.
    y=(2+PMLw):((SizeJ+2*PMLw+1)-PMLw-1);
    Bx(x,y,nf) = Bx(x,y,n0) + (-Ez(x,y,n0) + Ez(x,y-1,n0)) * dt/delta;
    Hx(x,y,nf) = am(x,y).*(Bx(x,y,nf)-2*Bx(x,y,n0)+Bx(x,y,np))+bm(x,y).*(Bx(x,y,nf)-Bx(x,y,np))+cm(x,y).*(2*Hx(x,y,n0)-Hx(x,y,np))+dm(x,y).*(2*Hx(x,y,n0)+Hx(x,y,np))+em(x,y).*Hx(x,y,np);
    if PMLw > 0
        % Bx in lower PML layer.
        y=2:PMLw+1;
        Bx(x,y,nf) = Bx(x,y,n0) + dt*((1/kappmy)*(-Ez(x,y,n0) + Ez(x,y-1,n0)) * 1/delta + PsiHxY(x,y));
        Hx(x,y,nf) = Bx(x,y,nf)./(u0*uinf(x,y));
        % Bx in upper PML layer.
        y=(SizeJ+2*PMLw+1)-PMLw:(SizeJ+2*PMLw);
        Bx(x,y,nf) = Bx(x,y,n0) + dt*((1/kappmy)*(-Ez(x,y,n0) + Ez(x,y-1,n0)) * 1/delta + PsiHxY(x,y));
        Hx(x,y,nf) = Bx(x,y,nf)./(u0*uinf(x,y));
    end

    % ========================= By and Hy =============================
    % Hy Psi array.
    x=1:SizeI-1;
    y=1:SizeJ+2*PMLw;
    PsiHyX(x,y) = (Cmx/delta)*(Ez(x+1,y,n0)-Ez(x,y,n0)) + bmx*PsiHyX(x,y);
    PsiHyX(SizeI,y) = (Cmx/delta)*(Ez(1,y,n0)-Ez(SizeI,y,n0)) + bmx*PsiHyX(SizeI,y);
    % By in normal space.
    y=(1+PMLw):(SizeJ+2*PMLw)-PMLw;
    By(x,y,nf) = By(x,y,n0) + (Ez(x+1,y,n0) - Ez(x,y,n0)) * dt/delta;
    By(SizeI,y,nf) = By(SizeI,y,n0) + (Ez(1,y,n0) - Ez(SizeI,y,n0)) * dt/delta; % PBC
    x=1:SizeI;
    Hy(x,y,nf) = am(x,y).*(By(x,y,nf)-2*By(x,y,n0)+By(x,y,np))+bm(x,y).*(By(x,y,nf)-By(x,y,np))+cm(x,y).*(2*Hy(x,y,n0)-Hy(x,y,np))+dm(x,y).*(2*Hy(x,y,n0)+Hy(x,y,np))+em(x,y).*Hy(x,y,np);
    if PMLw > 0
        % By in lower PML layer.
        x=1:SizeI-1;
        y=1:PMLw;
        By(x,y,nf) = By(x,y,n0) + dt*((1/kappmx)*(Ez(x+1,y,n0) - Ez(x,y,n0)) * 1/delta + PsiHyX(x,y));
        By(SizeI,y,nf) = By(SizeI,y,n0) + dt*((1/kappmx)*(Ez(1,y,n0) - Ez(SizeI,y,n0)) * 1/delta + PsiHyX(SizeI,y)); % PBC
        x=1:SizeI;
        Hy(x,y,nf) = By(x,y,nf)./(u0*uinf(x,y));
        % By in upper PML layer.
        x=1:SizeI-1;
        y=(SizeJ+2*PMLw)-PMLw+1:(SizeJ+2*PMLw);
        By(x,y,nf) = By(x,y,n0) + dt*((1/kappmx)*(Ez(x+1,y,n0) - Ez(x,y,n0)) * 1/delta + PsiHyX(x,y));
        By(SizeI,y,nf) = By(SizeI,y,n0) + dt*((1/kappmx)*(Ez(1,y,n0) - Ez(SizeI,y,n0)) * 1/delta + PsiHyX(SizeI,y)); % PBC
        x=1:SizeI;
        Hy(x,y,nf) = By(x,y,nf)./(u0*uinf(x,y));
    end

    % ========================= Dz and Ez =============================
    % Psi arrays.
    x=2:SizeI;
    y=1:SizeJ+2*PMLw;
    PsiEzX(x,y) = (Cex/delta)*(Hy(x,y,nf)-Hy(x-1,y,nf)) + bex*PsiEzX(x,y);
    PsiEzX(1,y) = (Cex/delta)*(Hy(1,y,nf)-Hy(SizeI,y,nf)) + bex*PsiEzX(1,y); % PBC
    PsiEzY(x,y) = (Cey/delta)*(-Hx(x,y+1,nf)+Hx(x,y,nf)) + bey*PsiEzY(x,y);
    PsiEzY(1,y) = (Cey/delta)*(-Hx(1,y+1,nf)+Hx(1,y,nf)) + bey*PsiEzY(1,y); % PBC
    % Dz in Normal Space.
    y=(1+PMLw):((SizeJ+2*PMLw)-PMLw);
    Dz(x,y,nf) = Dz(x,y,n0) + (Hy(x,y,nf)-Hy(x-1,y,nf)-Hx(x,y+1,nf)+Hx(x,y,nf)) * dt/delta;
    Dz(1,y,nf) = Dz(1,y,n0) + (Hy(1,y,nf)-Hy(SizeI,y,nf)-Hx(1,y+1,nf)+Hx(1,y,nf)) * dt/delta; % PBC
    x=1:SizeI;
    Ez(x,y,nf) = ae(x,y).*(Dz(x,y,nf)-2*Dz(x,y,n0)+Dz(x,y,np))+be(x,y).*(Dz(x,y,nf)-Dz(x,y,np))+ce(x,y).*(2*Ez(x,y,n0)-Ez(x,y,np))+de(x,y).*(2*Ez(x,y,n0)+Ez(x,y,np))+ee(x,y).*Ez(x,y,np);
    if PMLw > 0
        % Dz in lower PML layer.
        x=2:SizeI;
        y=1:PMLw;
        Dz(x,y,nf) = Dz(x,y,n0) + dt*(((1/kappex)*(Hy(x,y,nf)-Hy(x-1,y,nf))+(1/kappey)*(-Hx(x,y+1,nf)+Hx(x,y,nf))) * 1/delta + PsiEzX(x,y) + PsiEzY(x,y));
        Dz(1,y,nf) = Dz(1,y,n0) + dt*(((1/kappex)*(Hy(1,y,nf)-Hy(SizeI,y,nf))+(1/kappey)*(-Hx(1,y+1,nf)+Hx(1,y,nf))) * 1/delta + PsiEzX(1,y) + PsiEzY(1,y)); % PBC
        x=1:SizeI;
        Ez(x,y,nf) = Dz(x,y,nf)./(e0*einf(x,y));
        % Dz in upper PML layer.
        x=2:SizeI;
        y=(SizeJ+2*PMLw)-PMLw+1:(SizeJ+2*PMLw);
        Dz(x,y,nf) = Dz(x,y,n0) + dt*(((1/kappex)*(Hy(x,y,nf)-Hy(x-1,y,nf))+(1/kappey)*(-Hx(x,y+1,nf)+Hx(x,y,nf))) * 1/delta + PsiEzX(x,y) + PsiEzY(x,y));
        Dz(1,y,nf) = Dz(1,y,n0) + dt*(((1/kappex)*(Hy(1,y,nf)-Hy(SizeI,y,nf))+(1/kappey)*(-Hx(1,y+1,nf)+Hx(1,y,nf))) * 1/delta + PsiEzX(1,y) + PsiEzY(1,y)); % PBC
        x=1:SizeI;
        Ez(x,y,nf) = Dz(x,y,nf)./(e0*einf(x,y));
    end

    % ====================== Source ===================
    if SourcePlane == 1
        x = 1:SizeI;
        y = SourceLocationY;        
    else
        x = SourceLocationX;
        y = SourceLocationY;
    end
    % Source.
    if SourceChoice == 1
        Ez(x,y,nf) = Ez(x,y,nf) + exp( -1*((n-td)/(PulseWidth/4))^2 ) * Sc;
    elseif SourceChoice == 2
        Ez(x,y,nf) = Ez(x,y,nf) + sin(2*pi*f*(n)*dt) * Sc;
    elseif SourceChoice == 3
        Ez(x,y,nf) = Ez(x,y,nf) + (1-2*(pi*fp*(n*dt-dr))^2)*exp(-1*(pi*fp*(n*dt-dr))^2) * Sc;
    end
    Dz(x,y,nf) = e0*Ez(x,y,nf);

    % Transmitted fields.
    Ezt(n+1) = Ez(SizeI/2,SlabLeft+1,nf);
    Eztt(n+1) = Ez(SizeI/2,SlabRight+10,nf);

    % Fields for calculation of refractive index.
    Ezy1(n+1) = Ez(SizeI/2,Y1,nf);
    Ezy2(n+1) = Ez(SizeI/2,Y2, nf);
    
    if (mod(n, SnapshotInterval) == 0)
        EzSnapshots(:,:,n/SnapshotInterval+1) = Ez(1+(0:SnapshotResolution:(SizeI-1)), 1+(0:SnapshotResolution:((SizeJ+2*PMLw)-1)), nf);
    end

    np = mod(np, 3)+1;
    n0 = mod(n0, 3)+1;
    nf = mod(nf, 3)+1;
end
fprintf ( 1, '\rSimulation complete! \n');
toc
% Electric field snapshots.
for i=1:(MaxTime/SnapshotInterval)-1

    figure (6)
    mesh ( EzSnapshots (:, :, i) );
    view (4, 4)
    zlim ( [-1 1] )
    caxis([-0.1 0.6])
    xlabel ('y-axis')
    ylabel ('x-axis')
    %colorbar

    figure (7)
    mesh ( EzSnapshots (:, :, i) );
    view (0, 90)
    zlim ( [-10 10] )
    caxis([-0.1 0.6])
    xlabel ('y-axis')
    ylabel ('x-axis')
    %colorbar

end

% Postprocessing.
Fs = 1/dt;                    % Sampling frenuency
T = dt;                       % Sample time
L = length(Ezi);              % Length of signal
t = (0:L-1)*T;                % Time vector
fspan = 100;                  % Points to plot in frenuency domain

figure(1)
subplot(211)
plot(Fs*t, Ezi, 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Incident Field', 'FontSize', 12, 'FontWeight', 'b')
xlabel('time step (n)', 'FontSize', 11, 'FontWeight', 'b')
grid on
figure(2)
subplot(211)
plot(Fs*t, Ezt, 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Transmitted Field', 'FontSize', 12, 'FontWeight', 'b')
xlabel('time step (n)', 'FontSize', 11, 'FontWeight', 'b')
grid on
figure(3)
subplot(211)
plot(Fs*t, Eztt, 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Transmitted Field Beyond Slab', 'FontSize', 12, 'FontWeight', 'b')
xlabel('time step (n)', 'FontSize', 11, 'FontWeight', 'b')
grid on

NFFT = 2^nextpow2(L); % Next power of 2 from length of Exi
% Incident and Transmitted fields.
EZI = fft(Ezi,NFFT)/L;
EZT = fft(Ezt,NFFT)/L;
EZTT = fft(Eztt,NFFT)/L;
% Refractive index calculations.
EZY1 = fft(Ezy1,NFFT)/L;
EZY2 = fft(Ezy2,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure(1)
subplot(212)
EZIp = 2*abs(EZI(1:NFFT/2+1));
plot(f(1:fspan), EZIp(1:fspan), 'LineWidth', 2.0, 'Color', 'r')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Single-Sided Amplitude Spectrum of Ezi(t)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frenuency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('|EZI(f)|', 'FontSize', 11, 'FontWeight', 'b')
grid on
figure(2)
subplot(212)
EZTp = 2*abs(EZT(1:NFFT/2+1));
plot(f(1:fspan), EZTp(1:fspan), 'LineWidth', 2.0, 'Color', 'r')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Single-Sided Amplitude Spectrum of Ezt(t)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frenuency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('|EZT(f)|', 'FontSize', 11, 'FontWeight', 'b')
grid on
figure(3)
subplot(212)
EZTTp = 2*abs(EZTT(1:NFFT/2+1));
plot(f(1:fspan), EZTTp(1:fspan), 'LineWidth', 2.0, 'Color', 'r')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Single-Sided Amplitude Spectrum of Eztt(t)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frenuency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('|EZT(f)|', 'FontSize', 11, 'FontWeight', 'b')
grid on

% Transmission Coefficient.
figure(4)
subplot(211)
TAU = abs(EZT(1:NFFT/2+1)./EZI(1:NFFT/2+1));
plot(f(1:fspan), TAU(1:fspan), 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Transmission Coefficient', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frenuency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('|EZT(f)/EZI(f)|', 'FontSize', 11, 'FontWeight', 'b')
axis([-1 1 -2 2])
axis 'auto x'
grid on
subplot(212)
plot(f(1:fspan), 1-TAU(1:fspan), 'LineWidth', 2.0, 'Color', 'b')
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Reflection Coefficient', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frenuency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('1-|EZT(f)/EZI(f)|', 'FontSize', 11, 'FontWeight', 'b')
axis([-1 1 -2 2])
axis 'auto x'
grid on

% Refractive Index calculations.
nFDTD = (1/(1i*k0*(y1-y2))).*log(EZY2(1:NFFT/2+1)./EZY1(1:NFFT/2+1));
figure(5)
subplot(211)
plot(f(1:fspan), real(nFDTD(1:fspan)), 'LineWidth', 2.0, 'Color', 'b');
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Refractive index re(n)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frenuency (Hz)', 'FontSize', 11)
ylabel('re(n)', 'FontSize', 11)
grid on
subplot(212)
plot(f(1:fspan), imag(nFDTD(1:fspan)), 'LineWidth', 2.0, 'Color', 'r');
set(gca, 'FontSize', 10, 'FontWeight', 'b')
title('Refractive index im(n)', 'FontSize', 12, 'FontWeight', 'b')
xlabel('Frenuency (Hz)', 'FontSize', 11, 'FontWeight', 'b')
ylabel('im(n)', 'FontSize', 11, 'FontWeight', 'b')
grid on
\end{lstlisting}
\chapter{GPU Implementation of 1D and 2D DNG}
\section{1D DNG Class Definition for CUDA Implementation}
\label{App:1D-DNG-Class-Definition-CUDA}
\lstset{language=[ISO]C++, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black}}
\begin{lstlisting}
class CFDTD1DDNG
{
private:
	// Simulation parameters.
	const unsigned int Size;
	const unsigned int MaxTime;
	const unsigned int PulseWidth;
	const unsigned int td;
	const unsigned int SourceLocation;
	const unsigned int SlabLeft;
	const unsigned int SlabRight;
	const unsigned int SnapshotInterval;

	// Choice of source.
	// 1. Gaussian pulse 2. Sine wave 3. Ricker wavelet
	const unsigned int SourceChoice;

	const PRECISION e0;
	const PRECISION u0;
	const PRECISION dt;
	const PRECISION dz;
	const PRECISION Sc;

	// Frequency, wavelength, wave number.
	const PRECISION l;
	const PRECISION f;
	const PRECISION fmax;
	const PRECISION w;
	const PRECISION k0;
	const PRECISION fp; // Ricker wavelet peak frequency.
	const PRECISION dr; // Ricker wavelet delay.

	// Data arrays.
	PRECISION *Ex_;
	PRECISION *Dx_;
	PRECISION *Hy_;
	PRECISION *By_;
	unsigned int frame;

	// Incident and transmitted fields.
	PRECISION *Exi;
	PRECISION *Ext;
	PRECISION *Extt;
	const unsigned int x1;

	// Refractive index.
	const unsigned int Z1;
	const PRECISION z1;
	const unsigned int Z2;
	const PRECISION z2;
	PRECISION *Exz1;
	PRECISION *Exz2;

	// Drude material parameters.
	PRECISION *einf;
	PRECISION *uinf;
	PRECISION *wpesq;
	PRECISION *wpmsq;
	PRECISION *ge;
	PRECISION *gm;

	// Auxiliary field scalars.
	PRECISION *ae0, *ae, *be, *ce, *de, *ee;
	PRECISION *am0, *am, *bm, *cm, *dm, *em;

	// Time indices.
	unsigned int nf, n0, np;

	// Timer variables.
	__int64 tStart;
	__int64 tEnd;

	// ================== Device arrays ==================
	// Data arrays.
	PRECISION *d_Ex_;
	PRECISION *d_Dx_;
	PRECISION *d_Hy_;
	PRECISION *d_By_;
	// Incident and transmitted fields.
	PRECISION *d_Exi;
	PRECISION *d_Ext;
	PRECISION *d_Extt;
	// Refractive Index.
	PRECISION *d_Exz1;
	PRECISION *d_Exz2;
	// Drude material parameters.
	PRECISION *d_einf;
	PRECISION *d_uinf;
	PRECISION *d_wpesq;
	PRECISION *d_wpmsq;
	PRECISION *d_ge;
	PRECISION *d_gm;
	// Auxiliary field scalars.
	PRECISION *d_ae0, *d_ae, *d_be, *d_ce, *d_de, *d_ee;
	PRECISION *d_am0, *d_am, *d_bm, *d_cm, *d_dm, *d_em;
	// ===================================================

public:
	CFDTD1DDNG(unsigned int=4U*1024U, unsigned int=10U, unsigned int=16U, unsigned int=1U);

	// Space calculations.
	unsigned long SimSize();
	unsigned long HDDSpace();

	// Memory allocation and initialisation.
	int AllocateMemoryCPU();
	int InitialiseCPU();
	int InitialiseExHyCPU();
	int AllocateMemoryGPU();
	int CopyDataCPUtoGPU();
	int CopyExHyCPUtoGPU();

	// Simulations.
	int DryRunCPU();
	int RunSimulationCPU(bool=true);
	int DryRunGPU();
	int RunSimulationGPU(bool=true);

	// Complete Runs on CPU and GPU.
	int CompleteRunCPU(bool=true);
	int CompleteRunGPU(bool=true);

	// Timing.
	void StartTimer();
	void StopTimer();
	PRECISION GetElapsedTime();

	int SafeCall(int, const char[]=NULL);

	int CleanupCPU();
	int CleanupGPU();
	~CFDTD1DDNG();
};
\end{lstlisting}
\section{2D DNG Class Definition for OpenCL Implementation}
\label{App:2D-DNG-Class-Definition-OpenCL}
\begin{lstlisting}
class CFDTD2DDNG
{
private:
	// Simulation parameters.
	const unsigned int I;
	const unsigned int J;
	const unsigned int PMLw;
	const unsigned int SlabLeft;
	const unsigned int SlabRight;
	const unsigned int MaxTime;
	const unsigned int PulseWidth;
	const unsigned int td;
	const unsigned int SnapshotResolution;
	const unsigned int SnapshotInterval;
	// Choice of source.
	// 1. Gaussian pulse 2. Sine wave 3. Ricker wavelet
	const unsigned int SourceChoice;
	const unsigned int SourcePlane;	// 0. Omni, 1. Plane wave.
	const unsigned int SourceLocationX;
	const unsigned int SourceLocationY;

	const PRECISION c;
	const PRECISION pi;
	const PRECISION e0;
	const PRECISION u0;
	const PRECISION dt;
	const PRECISION delta;
	const PRECISION Sc;

	// Frequency, wavelength, wave number.
	const PRECISION l;
	const PRECISION f;
	const PRECISION fmax;
	const PRECISION w;
	const PRECISION k0;
	const PRECISION fp; // Ricker wavelet peak frequency.
	const PRECISION dr; // Ricker wavelet delay.

	// Data array sizes.
	const unsigned int IEz, JEz;
	const unsigned int IHx, JHx;
	const unsigned int IHy, JHy;

	// Data arrays.
	PRECISION *Ez_;
	PRECISION *Dz_;
	PRECISION *Hx_;
	PRECISION *Bx_;
	PRECISION *Hy_;
	PRECISION *By_;

	// Incident and transmitted fields.
	PRECISION *Ezi;
	PRECISION *Ezt;
	PRECISION *Eztt;
	const unsigned int x1;

	// Refractive index.
	const unsigned int Y1;
	const PRECISION y1;
	const unsigned int Y2;
	const PRECISION y2;
	PRECISION *Ezy1;
	PRECISION *Ezy2;

	// Drude material parameters.
	PRECISION *einf_;
	PRECISION *uinf_;
	PRECISION *wpesq_;
	PRECISION *wpmsq_;
	PRECISION *ge_;
	PRECISION *gm_;

	// Auxiliary field scalars.
	PRECISION *ae0_, *ae_, *be_, *ce_, *de_, *ee_;
	PRECISION *am0_, *am_, *bm_, *cm_, *dm_, *em_;

	// PML arrays.
	PRECISION *PsiEzX_, *PsiEzY_;
	PRECISION *PsiHyX_, *PsiHxY_;

	// PML parameters.
	const PRECISION kappe, kappm;
	const PRECISION kappex, kappey, kappmx, kappmy;
	const PRECISION aex, aey, amx, amy;
	const PRECISION sigex, sigey, sigmx, sigmy;
	const PRECISION bex, bey, bmx, bmy;
	const PRECISION Cex, Cey, Cmx, Cmy;

	// Snapshot frame number.
	unsigned int frame;

	// Time indices.
	unsigned int nf, n0, np;

	// Timer variables.
	__int64 tStart;
	__int64 tEnd;

	// ================== Device arrays ==================
	// Data arrays.
	cl_mem d_Ez_;
	cl_mem d_Dz_;
	cl_mem d_Hx_;
	cl_mem d_Bx_;
	cl_mem d_Hy_;
	cl_mem d_By_;
	// Incident and transmitted fields.
	cl_mem d_Ezi;
	cl_mem d_Ezt;
	cl_mem d_Eztt;
	// Refractive Index.
	cl_mem d_Ezy1;
	cl_mem d_Ezy2;
	// Drude material parameters.
	cl_mem d_einf_;
	cl_mem d_uinf_;
	cl_mem d_wpesq_;
	cl_mem d_wpmsq_;
	cl_mem d_ge_;
	cl_mem d_gm_;
	// Auxiliary field scalars.
	cl_mem d_ae0_, d_ae_, d_be_, d_ce_, d_de_, d_ee_;
	cl_mem d_am0_, d_am_, d_bm_, d_cm_, d_dm_, d_em_;
	// PML arrays.
	cl_mem d_PsiEzX_, d_PsiEzY_;
	cl_mem d_PsiHyX_, d_PsiHxY_;
	// ===================================================

	// OPENCL related parameters.
	bool cpu;		// Is OPENCL using CPU or GPU?

	// OPENCL context/device/program
	cl_context context;
	cl_device_id *devices;
	cl_command_queue commandQueue;
	cl_program program;

	// Dry run and simulation kernel handles.
	cl_kernel DryRun_kernel_M;
	cl_kernel DryRun_kernel_E;
	cl_kernel Simulation_kernel_M;
	cl_kernel Simulation_kernel_E;

public:
	CFDTD2DDNG(
				unsigned int=512,	// I
				unsigned int=512,	// J
				unsigned int=64,	// PMLw
				unsigned int=4*512,	// MaxTime
				unsigned int=1,		// Snapshot resolution
				unsigned int=10,	// Snapshot interval
				unsigned int=1,		// Source choice
				unsigned int=1,		// Source is plane wave?
				unsigned int=256,	// Source location X
				unsigned int=64+5);	// Source location Y

	// Space calculations.
	unsigned long SimSize();
	unsigned long HDDSpace();

	// Memory allocation and initialisation.
	int AllocateMemoryCPU();
	int InitialiseCPU();
	int InitialiseCL();				// Search and allocate a device.
	int InitialiseForSimulationCPU();
	int AllocateMemoryGPU();
	int InitialiseCLKernelsGPU();	// Build/attach kernels to respective kernel functions and set arguments.
	int InitialiseForSimulationGPU();

	// Simulations.
	int DryRunCPU();
	int RunSimulationCPU(bool=true);
	int DryRunGPU();
	int RunSimulationGPU(bool=true);

	// Complete Runs on CPU and GPU.
	int CompleteRunCPU(bool=true);
	int CompleteRunGPU(bool=true);

	// Timing.
	void StartTimer();
	void StopTimer();
	PRECISION GetElapsedTime();

	std::string convertToString(const char * filename);

	int SafeCall(int, const char []=NULL);

	int CleanupCPU();
	int CleanupCL();
	int CleanupGPU();
	~CFDTD2DDNG ();
};
\end{lstlisting}
\section{CUDA 1D DNG Kernels}
\label{App:CUDA-1D-DNG-Kernels}
\begin{lstlisting}
#define PRECISION double

// Simulation kernels.
template <unsigned int BlockX, unsigned int BlockY> __global__ void FDTD1DDNGKernel_Simulation_M(
							const unsigned int Size,
							const unsigned int PulseWidth,
							const unsigned int td,
							const unsigned int SourceLocation,
							const unsigned int SourceChoice,
							const PRECISION e0,
							const PRECISION u0,
							const PRECISION dt,
							const PRECISION dz,
							const PRECISION Sc,
							// Frequency, wavelength, wave number.
							const PRECISION f,
							const PRECISION fp,
							const PRECISION dr,
							// Data arrays.
							PRECISION *Ex_, PRECISION *Dx_, PRECISION *Hy_, PRECISION *By_,
							// Drude material parameters.
							PRECISION *einf, PRECISION *uinf, PRECISION *wpesq, PRECISION *wpmsq, PRECISION *ge, PRECISION *gm,
							// Drude scalars.
							PRECISION *ae0, PRECISION *ae, PRECISION *be, PRECISION *ce, PRECISION *de, PRECISION *ee,
							PRECISION *am0, PRECISION *am, PRECISION *bm, PRECISION *cm, PRECISION *dm, PRECISION *em,
							// Incident field.
							PRECISION *Ext,
							PRECISION *Extt,
							PRECISION *Exz1,
							PRECISION *Exz2,
							const unsigned int x1,
							const unsigned int Z1,
							const unsigned int Z2,
							// Time indices.
							const unsigned int n,
							const unsigned int np,
							const unsigned int n0,
							const unsigned int nf)
{
	unsigned int i = BlockX*blockIdx.x+threadIdx.x;

	if (i != Size-1) // Normal update equation.
	{
		By(i,nf) = By(i,n0) + (Ex(i,n0)-Ex(i+1,n0))*dt/dz;
		Hy(i,nf) = am[i]*(By(i,nf)-2*By(i,n0)+By(i,np)) + bm[i]*(By(i,nf)-By(i,np)) + cm[i]*(2*Hy(i,n0)-Hy(i,np)) + dm[i]*(2*Hy(i,n0)+Hy(i,np)) + em[i]*(Hy(i,np));
	}
	__syncthreads();

	// ABC
	if (i == Size-1)
	{
		Hy(i,nf) = Hy(i-1,n0) + (Sc-1)/(Sc+1)*(Hy(i-1,nf)-Hy(i,n0));
		By(i,nf) = u0*Hy(i,nf);
	}
}
template <unsigned int BlockX, unsigned int BlockY> __global__ void FDTD1DDNGKernel_Simulation_E(
							const unsigned int Size,
							const unsigned int PulseWidth,
							const unsigned int td,
							const unsigned int SourceLocation,
							const unsigned int SourceChoice,
							const PRECISION e0,
							const PRECISION u0,
							const PRECISION dt,
							const PRECISION dz,
							const PRECISION Sc,
							// Frequency, wavelength, wave number.
							const PRECISION f,
							const PRECISION fp,
							const PRECISION dr,
							// Data arrays.
							PRECISION *Ex_, PRECISION *Dx_, PRECISION *Hy_, PRECISION *By_,
							// Drude material parameters.
							PRECISION *einf, PRECISION *uinf, PRECISION *wpesq, PRECISION *wpmsq, PRECISION *ge, PRECISION *gm,
							// Drude scalars.
							PRECISION *ae0, PRECISION *ae, PRECISION *be, PRECISION *ce, PRECISION *de, PRECISION *ee,
							PRECISION *am0, PRECISION *am, PRECISION *bm, PRECISION *cm, PRECISION *dm, PRECISION *em,
							// Incident field.
							PRECISION *Ext,
							PRECISION *Extt,
							PRECISION *Exz1,
							PRECISION *Exz2,
							const unsigned int x1,
							const unsigned int Z1,
							const unsigned int Z2,
							// Time indices.
							const unsigned int n,
							const unsigned int np,
							const unsigned int n0,
							const unsigned int nf)
{
	unsigned int i = BlockX*blockIdx.x+threadIdx.x;

	if (i != 0)
	{
		Dx(i,nf) = Dx(i,n0) + (Hy(i-1,nf)-Hy(i,nf))*dt/dz;
		Ex(i,nf) = ae[i]*(Dx(i,nf)-2*Dx(i,n0)+Dx(i,np)) + be[i]*(Dx(i,nf)-Dx(i,np)) + ce[i]*(2*Ex(i,n0)-Ex(i,np)) + de[i]*(2*Ex(i,n0)+Ex(i,np)) + ee[i]*(Ex(i,np));
	}
	__syncthreads();

	// ABC
	if (i == 0)
	{
		Ex(i,nf) = Ex(i+1,n0) + (Sc-1)/(Sc+1)*(Ex(i+1,nf)-Ex(i,n0));
		Dx(i,nf) = e0*Ex(i,nf);
	}
	__syncthreads();

	// Source.
	if (i == SourceLocation)
	{
		if (SourceChoice == 1)
		{
			Ex(SourceLocation,nf) = Ex(SourceLocation,nf) + exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
		}
		else if (SourceChoice == 2)
		{
			Ex(SourceLocation,nf) = Ex(SourceLocation,nf) + sin(2.*PI*f*(PRECISION)n*dt) * Sc;
		}
		else if (SourceChoice == 3)
		{
			Ex(SourceLocation,nf) = Ex(SourceLocation,nf) + (1.-2.*pow(PI*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(PI*fp*((PRECISION)n*dt-dr),2)) * Sc;
		}
		Dx(SourceLocation,nf) = e0*Ex(SourceLocation,nf);
	}
	// Recording transmitted Fields.
	if (i == x1)
		Ext[n] = Ex(i,nf);
	if (i == (Size-(2*Size/3))+10)
		Extt[n] = Ex(i,nf);
	if (i == Z1)
		Exz1[n] = Ex(i,nf);
	if (i == Z2)
		Exz2[n] = Ex(i,nf);
}
\end{lstlisting}
\section{OpenCL 2D DNG Kernels}
\label{App:OpenCL-2D-DNG-Kernels}
\begin{lstlisting}
#define PRECISION double

#ifdef cl_amd_fp64
#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#else
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#endif

// Number of const uint or double arguments does NOT have any significant impact on performance.
__kernel void FDTD2DDNGKernel_DryRun_M(
							const unsigned int I,
							const unsigned int J,
							const unsigned int PMLw,
							const unsigned int PulseWidth,
							const unsigned int td,
							const unsigned int SourceChoice,
							const unsigned int SourcePlane,
							const unsigned int SourceLocationX,
							const unsigned int SourceLocationY,
							const PRECISION c,
							const PRECISION pi,
							const PRECISION e0,
							const PRECISION u0,
							const PRECISION dt,
							const PRECISION delta,
							const PRECISION Sc,
							const PRECISION f,
							const PRECISION fp,
							const PRECISION dr,
							const unsigned int IEz, const unsigned int JEz,
							const unsigned int IHx, const unsigned int JHx,
							const unsigned int IHy, const unsigned int JHy,
							const __global PRECISION *Ez_, const __global PRECISION *Dz_, __global PRECISION *Hx_, __global PRECISION *Bx_, __global PRECISION *Hy_, __global PRECISION *By_,
							const __global PRECISION *PsiEzX_, const __global PRECISION *PsiEzY_, __global PRECISION *PsiHyX_, __global PRECISION *PsiHxY_,
							const __global PRECISION *einf_, const __global PRECISION *uinf_,
							const PRECISION kappex, const PRECISION kappey, const PRECISION kappmx, const PRECISION kappmy,
							const PRECISION bex, const PRECISION bey, const PRECISION bmx, const PRECISION bmy,
							const PRECISION Cex, const PRECISION Cey, const PRECISION Cmx, const PRECISION Cmy,
							__global PRECISION *Ezi,
							const unsigned int x1,
							const unsigned int n,
							const unsigned int np,
							const unsigned int n0,
							const unsigned int nf)
{
	const unsigned int i = get_global_id(0);
	const unsigned int j = get_global_id(1);

	// Calculation of PsiHxY.
	if (j>0 && j<JHx-1)
		PsiHxY(i,j) = (Cmy/delta)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) + bmy*PsiHxY(i,j);

	// Bx in normal space.
	if( j>PMLw && j<JHx-PMLw-1)
	{
		Bx(i,j,nf) = Bx(i,j,n0) + (-Ez(i,j,n0) + Ez(i,j-1,n0)) * dt/delta;
		Hx(i,j,nf) = Bx(i,j,nf)/(u0*uinf(i,j));
	}
	// Bx in lower PML.
	if (j>0 && j<PMLw+1)
	{
		Bx(i,j,nf) = Bx(i,j,n0) + dt*((1./kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1./delta + PsiHxY(i,j));
		Hx(i,j,nf) = Bx(i,j,nf)/(u0*uinf(i,j));
	}
	// Bx in upper PML.
	if (j>JHx-PMLw-2 && j<JHx-1)
	{
		Bx(i,j,nf) = Bx(i,j,n0) + dt*((1./kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1./delta + PsiHxY(i,j));
		Hx(i,j,nf) = Bx(i,j,nf)/(u0*uinf(i,j));
	}

	// ========================== By and Hy ==========================
	if (i<IHy-1)
	{
		// PsiHyX array.
		if (j<JHy)
			PsiHyX(i,j) = (Cmx/delta)*(Ez(i+1,j,n0)-Ez(i,j,n0)) + bmx*PsiHyX(i,j);

		// By in normal space.
		if (j>PMLw-1 && j<JHy-PMLw)
		{
			By(i,j,nf) = By(i,j,n0) + (Ez(i+1,j,n0) - Ez(i,j,n0)) * dt/delta;
			Hy(i,j,nf) = By(i,j,nf)/(u0*uinf(i,j));
		}
		// By in Lower PML.
		if (j<PMLw)
		{
			By(i,j,nf) = By(i,j,n0) + dt*((1./kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1./delta + PsiHyX(i,j));
			Hy(i,j,nf) = By(i,j,nf)/(u0*uinf(i,j));
		}
		// By in upper PML.
		if (j>JHy-PMLw-1 && j<JHy)
		{
			By(i,j,nf) = By(i,j,n0) + dt*((1./kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1./delta + PsiHyX(i,j));
			Hy(i,j,nf) = By(i,j,nf)/(u0*uinf(i,j));
		}
	}
	else
	{
		// PsiHyX array.
		if (j<JHy)
			PsiHyX(IHy-1,j) = (Cmx/delta)*(Ez(0,j,n0)-Ez(IHy-1,j,n0)) + bmx*PsiHyX(IHy-1,j); // PBC

		// By in normal space.
		if (j>PMLw-1 && j<JHy-PMLw)
		{
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + (Ez(0,j,n0) - Ez(IHy-1,j,n0)) * dt/delta; // PBC
			Hy(IHy-1,j,nf) = By(IHy-1,j,nf)/(u0*uinf(IHy-1,j)); // PBC
		}
		// By in Lower PML.
		if (j<PMLw)
		{
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + dt*((1./kappmx)*(Ez(0,j,n0) - Ez(IHy-1,j,n0)) * 1./delta + PsiHyX(IHy-1,j)); // PBC
			Hy(IHy-1,j,nf) = By(IHy-1,j,nf)/(u0*uinf(IHy-1,j)); // PBC
		}
		// By in upper PML.
		if (j>JHy-PMLw-1 && j<JHy)
		{
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + dt*((1./kappmx)*(Ez(0,j,n0) - Ez(IHy-1,j,n0)) * 1./delta + PsiHyX(IHy-1,j)); // PBC
			Hy(IHy-1,j,nf) = By(IHy-1,j,nf)/(u0*uinf(IHy-1,j)); // PBC
		}
	}
}
__kernel void FDTD2DDNGKernel_DryRun_E(
							const unsigned int I,
							const unsigned int J,
							const unsigned int PMLw,
							const unsigned int PulseWidth,
							const unsigned int td,
							const unsigned int SourceChoice,
							const unsigned int SourcePlane,
							const unsigned int SourceLocationX,
							const unsigned int SourceLocationY,
							const PRECISION c,
							const PRECISION pi,
							const PRECISION e0,
							const PRECISION u0,
							const PRECISION dt,
							const PRECISION delta,
							const PRECISION Sc,
							const PRECISION f,
							const PRECISION fp,
							const PRECISION dr,
							const unsigned int IEz, const unsigned int JEz,
							const unsigned int IHx, const unsigned int JHx,
							const unsigned int IHy, const unsigned int JHy,
							__global PRECISION *Ez_, __global PRECISION *Dz_, const __global PRECISION *Hx_, const __global PRECISION *Bx_, const __global PRECISION *Hy_, const __global PRECISION *By_,
							__global PRECISION *PsiEzX_, __global PRECISION *PsiEzY_, const __global PRECISION *PsiHyX_, const __global PRECISION *PsiHxY_,
							const __global PRECISION *einf_, const __global PRECISION *uinf_,
							const PRECISION kappex, const PRECISION kappey, const PRECISION kappmx, const PRECISION kappmy,
							const PRECISION bex, const PRECISION bey, const PRECISION bmx, const PRECISION bmy,
							const PRECISION Cex, const PRECISION Cey, const PRECISION Cmx, const PRECISION Cmy,
							__global PRECISION *Ezi,
							const unsigned int x1,
							const unsigned int n,
							const unsigned int np,
							const unsigned int n0,
							const unsigned int nf)
{
	const unsigned int i = get_global_id(0);
	const unsigned int j = get_global_id(1);

	// ========================== Dz and Ez ==========================
	if (i>0)
	{
		// Psi arrays.
		if (j<JEz)
		{
			PsiEzX(i,j) = (Cex/delta)*(Hy(i,j,nf)-Hy(i-1,j,nf)) + bex*PsiEzX(i,j);
			PsiEzY(i,j) = (Cey/delta)*(-Hx(i,j+1,nf)+Hx(i,j,nf)) + bey*PsiEzY(i,j);
		}
		// Dz in normal space.
		if (j>PMLw-1 && j<JEz-PMLw)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + (Hy(i,j,nf)-Hy(i-1,j,nf)-Hx(i,j+1,nf)+Hx(i,j,nf)) * dt/delta;
			Ez(i,j,nf) = Dz(i,j,nf)/(e0*einf(i,j));
		}
		// Dz in lower PML.
		if (j<PMLw)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1./kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1./kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1./delta + PsiEzX(i,j) + PsiEzY(i,j));
			Ez(i,j,nf) = Dz(i,j,nf)/(e0*einf(i,j));
		}
		// Dz in upper PML.
		if (j>JEz-PMLw-1 && j<JEz)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1./kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1./kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1./delta + PsiEzX(i,j) + PsiEzY(i,j));
			Ez(i,j,nf) = Dz(i,j,nf)/(e0*einf(i,j));
		}
	}
	else
	{
		// Psi arrays.
		if (j<JEz)
		{
			PsiEzX(0,j) = (Cex/delta)*(Hy(0,j,nf)-Hy(IEz-1,j,nf)) + bex*PsiEzX(0,j); // PBC
			PsiEzY(0,j) = (Cey/delta)*(-Hx(0,j+1,nf)+Hx(0,j,nf)) + bey*PsiEzY(0,j); // PBC
		}
		// Dz in normal space.
		if (j>PMLw-1 && j<JEz-PMLw)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + (Hy(0,j,nf)-Hy(IEz-1,j,nf)-Hx(0,j+1,nf)+Hx(0,j,nf)) * dt/delta; // PBC
			Ez(0,j,nf) = Dz(0,j,nf)/(e0*einf(0,j)); // PBC
		}
		// Dz in lower PML.
		if (j<PMLw)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + dt*(((1./kappex)*(Hy(0,j,nf)-Hy(IEz-1,j,nf))+(1./kappey)*(-Hx(0,j+1,nf)+Hx(0,j,nf))) * 1./delta + PsiEzX(0,j) + PsiEzY(0,j)); // PBC
			Ez(0,j,nf) = Dz(0,j,nf)/(e0*einf(0,j)); // PBC
		}
		// Dz in upper PML.
		if (j>JEz-PMLw-1 && j<JEz)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + dt*(((1./kappex)*(Hy(0,j,nf)-Hy(IEz-1,j,nf))+(1./kappey)*(-Hx(0,j+1,nf)+Hx(0,j,nf))) * 1./delta + PsiEzX(0,j) + PsiEzY(0,j)); // PBC
			Ez(0,j,nf) = Dz(0,j,nf)/(e0*einf(0,j)); // PBC
		}
	}

	// ====================== Source ======================
	if (SourcePlane == 1 && i<IEz && j == SourceLocationY)
	{
		if (SourceChoice == 1)
			Ez(i,j,nf) += exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
		else if (SourceChoice == 2)
			Ez(i,j,nf) += sin(2.*pi*f*(PRECISION)n*dt) * Sc;
		else if (SourceChoice == 3)
			Ez(i,j,nf) += (1.-2.*pow(pi*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(pi*fp*((PRECISION)n*dt-dr),2)) * Sc;

		Dz(i,j,nf) = e0*Ez(i,j,nf);
	}
	else if (i == SourceLocationX && j == SourceLocationY)
	{
		if (SourceChoice == 1)
			Ez(i,j,nf) += exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
		else if (SourceChoice == 2)
			Ez(i,j,nf) += sin(2.*pi*f*(PRECISION)n*dt) * Sc;
		else if (SourceChoice == 3)
			Ez(i,j,nf) += (1.-2.*pow(pi*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(pi*fp*((PRECISION)n*dt-dr),2)) * Sc;

		Dz(i,j,nf) = e0*Ez(i,j,nf);
	}
	if (j==x1)
		Ezi[n] = Ez(IEz/2,x1,nf); // Incident field.
}
// Simulation kernel
__kernel void FDTD2DDNGKernel_Simulation_M(
							const unsigned int I,
							const unsigned int J,
							const unsigned int PMLw,
							const unsigned int PulseWidth,
							const unsigned int td,
							const unsigned int SourceChoice,
							const unsigned int SourcePlane,
							const unsigned int SourceLocationX,
							const unsigned int SourceLocationY,
							const PRECISION c,
							const PRECISION pi,
							const PRECISION e0,
							const PRECISION u0,
							const PRECISION dt,
							const PRECISION delta,
							const PRECISION Sc,
							const PRECISION f,
							const PRECISION fp,
							const PRECISION dr,
							const unsigned int IEz, const unsigned int JEz,
							const unsigned int IHx, const unsigned int JHx,
							const unsigned int IHy, const unsigned int JHy,
							const __global PRECISION *Ez_, const __global PRECISION *Dz_, __global PRECISION *Hx_, __global PRECISION *Bx_, __global PRECISION *Hy_, __global PRECISION *By_,
							const __global PRECISION *einf_, const __global PRECISION *uinf_, const __global PRECISION *wpesq_, const __global PRECISION *wpmsq_, const __global PRECISION *ge_, const __global PRECISION *gm_,
							const __global PRECISION *ae0_, const __global PRECISION *ae_, const __global PRECISION *be_, const __global PRECISION *ce_, const __global PRECISION *de_, const __global PRECISION *ee_,
							const __global PRECISION *am0_, const __global PRECISION *am_, const __global PRECISION *bm_, const __global PRECISION *cm_, const __global PRECISION *dm_, const __global PRECISION *em_,
							const __global PRECISION *PsiEzX_, const __global PRECISION *PsiEzY_, __global PRECISION *PsiHyX_, __global PRECISION *PsiHxY_,
							const PRECISION kappex, const PRECISION kappey, const PRECISION kappmx, const PRECISION kappmy,
							const PRECISION bex, const PRECISION bey, const PRECISION bmx, const PRECISION bmy,
							const PRECISION Cex, const PRECISION Cey, const PRECISION Cmx, const PRECISION Cmy,
							__global PRECISION *Ezt, __global PRECISION *Eztt, __global PRECISION *Ezy1, __global PRECISION *Ezy2,
							const unsigned int x1, const unsigned int Y1, const unsigned int Y2,
							const unsigned int n,
							const unsigned int np,
							const unsigned int n0,
							const unsigned int nf)
{
	const unsigned int i = get_global_id(0);
	const unsigned int j = get_global_id(1);

	// Calculation of PsiHxY.
	if (j>0 && j<JHx-1)
		PsiHxY(i,j) = (Cmy/delta)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) + bmy*PsiHxY(i,j);

	// Bx in normal space.
	if (j>PMLw && j<JHx-PMLw-1)
	{
		Bx(i,j,nf) = Bx(i,j,n0) + (-Ez(i,j,n0) + Ez(i,j-1,n0)) * dt/delta;
		Hx(i,j,nf) = am(i,j)*(Bx(i,j,nf)-2.*Bx(i,j,n0)+Bx(i,j,np))+bm(i,j)*(Bx(i,j,nf)-Bx(i,j,np))+cm(i,j)*(2.*Hx(i,j,n0)-Hx(i,j,np))+dm(i,j)*(2.*Hx(i,j,n0)+Hx(i,j,np))+em(i,j)*Hx(i,j,np);
	}
	// Bx in lower PML.
	if (j>0 && j<PMLw+1)
	{
		Bx(i,j,nf) = Bx(i,j,n0) + dt*((1./kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1./delta + PsiHxY(i,j));
		Hx(i,j,nf) = Bx(i,j,nf)/(u0*uinf(i,j));
	}
	// Bx in upper PML.
	if (j>JHx-PMLw-2 && j<JHx-1)
	{
		Bx(i,j,nf) = Bx(i,j,n0) + dt*((1./kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1./delta + PsiHxY(i,j));
		Hx(i,j,nf) = Bx(i,j,nf)/(u0*uinf(i,j));
	}


	// ========================== By and Hy ==========================
	if (i<IHy-1)
	{
		// PsiHyX arrays.
		if (j<JHy)
			PsiHyX(i,j) = (Cmx/delta)*(Ez(i+1,j,n0)-Ez(i,j,n0)) + bmx*PsiHyX(i,j);

		// By in normal space.
		if (j>PMLw-1 && j<JHy-PMLw)
		{
			By(i,j,nf) = By(i,j,n0) + (Ez(i+1,j,n0) - Ez(i,j,n0)) * dt/delta;
			Hy(i,j,nf) = am(i,j)*(By(i,j,nf)-2.*By(i,j,n0)+By(i,j,np))+bm(i,j)*(By(i,j,nf)-By(i,j,np))+cm(i,j)*(2.*Hy(i,j,n0)-Hy(i,j,np))+dm(i,j)*(2.*Hy(i,j,n0)+Hy(i,j,np))+em(i,j)*Hy(i,j,np);
		}
		// By in Lower PML.
		if (j<PMLw)
		{
			By(i,j,nf) = By(i,j,n0) + dt*((1./kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1./delta + PsiHyX(i,j));
			Hy(i,j,nf) = By(i,j,nf)/(u0*uinf(i,j));
		}
		// By in upper PML.
		if (j>JHy-PMLw-1 && j<JHy)
		{
			By(i,j,nf) = By(i,j,n0) + dt*((1./kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1./delta + PsiHyX(i,j));
			Hy(i,j,nf) = By(i,j,nf)/(u0*uinf(i,j));
		}
	}
	else
	{
		// PsiHyX array.
		if (j<JHy)
			PsiHyX(IHy-1,j) = (Cmx/delta)*(Ez(0,j,n0)-Ez(IHy-1,j,n0)) + bmx*PsiHyX(IHy-1,j); // PBC

		// By in normal space.
		if (j>PMLw-1 && j<JHy-PMLw)
		{
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + (Ez(0,j,n0) - Ez(IHy-1,j,n0)) * dt/delta; // PBC
			Hy(IHy-1,j,nf) = am(IHy-1,j)*(By(IHy-1,j,nf)-2.*By(IHy-1,j,n0)+By(IHy-1,j,np))+bm(IHy-1,j)*(By(IHy-1,j,nf)-By(IHy-1,j,np))+cm(IHy-1,j)*(2.*Hy(IHy-1,j,n0)-Hy(IHy-1,j,np))+dm(IHy-1,j)*(2.*Hy(IHy-1,j,n0)+Hy(IHy-1,j,np))+em(IHy-1,j)*Hy(IHy-1,j,np); // PBC
		}
		// By in Lower PML.
		if (j<PMLw)
		{
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + dt*((1./kappmx)*(Ez(0,j,n0) - Ez(IHy-1,j,n0)) * 1./delta + PsiHyX(IHy-1,j)); // PBC
			Hy(IHy-1,j,nf) = By(IHy-1,j,nf)/(u0*uinf(IHy-1,j)); // PBC
		}
		// By in upper PML.
		if (j>JHy-PMLw-1 && j<JHy)
		{
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + dt*((1./kappmx)*(Ez(0,j,n0) - Ez(IHy-1,j,n0)) * 1./delta + PsiHyX(IHy-1,j)); // PBC
			Hy(IHy-1,j,nf) = By(IHy-1,j,nf)/(u0*uinf(IHy-1,j)); // PBC
		}
	}
}
__kernel void FDTD2DDNGKernel_Simulation_E(
							const unsigned int I,
							const unsigned int J,
							const unsigned int PMLw,
							const unsigned int PulseWidth,
							const unsigned int td,
							const unsigned int SourceChoice,
							const unsigned int SourcePlane,
							const unsigned int SourceLocationX,
							const unsigned int SourceLocationY,
							const PRECISION c,
							const PRECISION pi,
							const PRECISION e0,
							const PRECISION u0,
							const PRECISION dt,
							const PRECISION delta,
							const PRECISION Sc,
							const PRECISION f,
							const PRECISION fp,
							const PRECISION dr,
							const unsigned int IEz, const unsigned int JEz,
							const unsigned int IHx, const unsigned int JHx,
							const unsigned int IHy, const unsigned int JHy,
							__global PRECISION *Ez_, __global PRECISION *Dz_, const __global PRECISION *Hx_, const __global PRECISION *Bx_, const __global PRECISION *Hy_, const __global PRECISION *By_,
							const __global PRECISION *einf_, const __global PRECISION *uinf_, const __global PRECISION *wpesq_, const __global PRECISION *wpmsq_, const __global PRECISION *ge_, const __global PRECISION *gm_,
							const __global PRECISION *ae0_, const __global PRECISION *ae_, const __global PRECISION *be_, const __global PRECISION *ce_, const __global PRECISION *de_, const __global PRECISION *ee_,
							const __global PRECISION *am0_, const __global PRECISION *am_, const __global PRECISION *bm_, const __global PRECISION *cm_, const __global PRECISION *dm_, const __global PRECISION *em_,
							__global PRECISION *PsiEzX_, __global PRECISION *PsiEzY_, const __global PRECISION *PsiHyX_, const __global PRECISION *PsiHxY_,
							const PRECISION kappex, const PRECISION kappey, const PRECISION kappmx, const PRECISION kappmy,
							const PRECISION bex, const PRECISION bey, const PRECISION bmx, const PRECISION bmy,
							const PRECISION Cex, const PRECISION Cey, const PRECISION Cmx, const PRECISION Cmy,
							__global PRECISION *Ezt, __global PRECISION *Eztt, __global PRECISION *Ezy1, __global PRECISION *Ezy2,
							const unsigned int x1, const unsigned int Y1, const unsigned int Y2,
							const unsigned int n,
							const unsigned int np,
							const unsigned int n0,
							const unsigned int nf)
{
	const unsigned int i = get_global_id(0);
	const unsigned int j = get_global_id(1);

	// ========================== Dz and Ez ==========================
	if (i>0)
	{
		// Psi arrays.
		if (j<JEz)
		{
			PsiEzX(i,j) = (Cex/delta)*(Hy(i,j,nf)-Hy(i-1,j,nf)) + bex*PsiEzX(i,j);
			PsiEzY(i,j) = (Cey/delta)*(-Hx(i,j+1,nf)+Hx(i,j,nf)) + bey*PsiEzY(i,j);
		}
		// Dz in normal space.
		if (j>PMLw-1 && j<JEz-PMLw)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + (Hy(i,j,nf)-Hy(i-1,j,nf)-Hx(i,j+1,nf)+Hx(i,j,nf)) * dt/delta;
			Ez(i,j,nf) = ae(i,j)*(Dz(i,j,nf)-2.*Dz(i,j,n0)+Dz(i,j,np))+be(i,j)*(Dz(i,j,nf)-Dz(i,j,np))+ce(i,j)*(2.*Ez(i,j,n0)-Ez(i,j,np))+de(i,j)*(2.*Ez(i,j,n0)+Ez(i,j,np))+ee(i,j)*Ez(i,j,np);
		}
		// Dz in lower PML.
		if (j<PMLw)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1./kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1./kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1./delta + PsiEzX(i,j) + PsiEzY(i,j));
			Ez(i,j,nf) = Dz(i,j,nf)/(e0*einf(i,j));
		}
		// Dz in upper PML.
		if (j>JEz-PMLw-1 && j<JEz)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1./kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1./kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1./delta + PsiEzX(i,j) + PsiEzY(i,j));
			Ez(i,j,nf) = Dz(i,j,nf)/(e0*einf(i,j));
		}
	}
	else
	{
		// Psi arrays.
		if (j<JEz)
		{
			PsiEzX(0,j) = (Cex/delta)*(Hy(0,j,nf)-Hy(IEz-1,j,nf)) + bex*PsiEzX(0,j); // PBC
			PsiEzY(0,j) = (Cey/delta)*(-Hx(0,j+1,nf)+Hx(0,j,nf)) + bey*PsiEzY(0,j); // PBC
		}
		// Dz in normal space.
		if (j>PMLw-1 && j<JEz-PMLw)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + (Hy(0,j,nf)-Hy(IEz-1,j,nf)-Hx(0,j+1,nf)+Hx(0,j,nf)) * dt/delta; // PBC
			Ez(0,j,nf) = ae(0,j)*(Dz(0,j,nf)-2.*Dz(0,j,n0)+Dz(0,j,np))+be(0,j)*(Dz(0,j,nf)-Dz(0,j,np))+ce(0,j)*(2.*Ez(0,j,n0)-Ez(0,j,np))+de(0,j)*(2.*Ez(0,j,n0)+Ez(0,j,np))+ee(0,j)*Ez(0,j,np); // PBC
		}
		// Dz in lower PML.
		if (j<PMLw)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + dt*(((1./kappex)*(Hy(0,j,nf)-Hy(IEz-1,j,nf))+(1./kappey)*(-Hx(0,j+1,nf)+Hx(0,j,nf))) * 1./delta + PsiEzX(0,j) + PsiEzY(0,j)); // PBC
			Ez(0,j,nf) = Dz(0,j,nf)/(e0*einf(0,j)); // PBC
		}
		// Dz in upper PML.
		if (j>JEz-PMLw-1 && j<JEz)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + dt*(((1./kappex)*(Hy(0,j,nf)-Hy(IEz-1,j,nf))+(1./kappey)*(-Hx(0,j+1,nf)+Hx(0,j,nf))) * 1./delta + PsiEzX(0,j) + PsiEzY(0,j)); // PBC
			Ez(0,j,nf) = Dz(0,j,nf)/(e0*einf(0,j)); // PBC
		}
	}

	// ====================== Source ======================
	if (SourcePlane == 1 && i<IEz && j == SourceLocationY)
	{
		if (SourceChoice == 1)
			Ez(i,j,nf) += exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
		else if (SourceChoice == 2)
			Ez(i,j,nf) += sin(2.*pi*f*(PRECISION)n*dt) * Sc;
		else if (SourceChoice == 3)
			Ez(i,j,nf) += (1.-2.*pow(pi*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(pi*fp*((PRECISION)n*dt-dr),2)) * Sc;

		Dz(i,j,nf) = e0*Ez(i,j,nf);
	}
	else if (i == SourceLocationX && j == SourceLocationY)
	{
		if (SourceChoice == 1)
			Ez(i,j,nf) += exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
		else if (SourceChoice == 2)
			Ez(i,j,nf) += sin(2.*pi*f*(PRECISION)n*dt) * Sc;
		else if (SourceChoice == 3)
			Ez(i,j,nf) += (1.-2.*pow(pi*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(pi*fp*((PRECISION)n*dt-dr),2)) * Sc;

		Dz(i,j,nf) = e0*Ez(i,j,nf);
	}
	// Recording transmitted fields.
	if (i==IEz/2 && j==x1)
		Ezt[n] = Ez(i,j,nf);
	if (i==IEz/2 && j==2*J/3+PMLw+1)
		Eztt[n] = Ez(i,j,nf);
	// Fields for refractive index.
	if (i==IEz/2 && j==Y1)
		Ezy1[n] = Ez(i,j,nf);
	if (i==IEz/2 && j==Y2)
		Ezy2[n] = Ez(i,j,nf);
}
\end{lstlisting}
\chapter{Cylindrical Cloak Implementations}
\section{Matlab Simulation Code}
\label{App:Matlab-Simulation-Cloak}
\lstset{language=Octave, commentstyle=\color{gray}, stringstyle=\color{red!60!black}}
\begin{lstlisting}
for n = 0:MaxTime
    
    % Progress indicator.
    if mod(n,SnapshotInterval) == 0
        fprintf(1, repmat('\b',1,linecount));
        linecount = fprintf(1, '%g %%', (n*100)/MaxTime );
    end
    
    % ========================= Bx =============================
    % Hx Psi array.
    i=1:IHx;
    j=2:JHx-1;
    PsiHxY(i,j) = (Cmy/delta)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) + bmy*PsiHxY(i,j);
    % Bx in normal space.
    j=(2+PMLw):(JHx-PMLw-1);
    Bx(i,j,nf) = Bx(i,j,n0) + (-Ez(i,j,n0) + Ez(i,j-1,n0)) * dt/delta;
    if PMLw > 0
        % Bx in lower PML layer.
        j=2:PMLw+1;
        Bx(i,j,nf) = Bx(i,j,n0) + dt*((1/kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1/delta + PsiHxY(i,j));
        % Bx in upper PML layer.
        j=JHx-PMLw:JHx-1;
        Bx(i,j,nf) = Bx(i,j,n0) + dt*((1/kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1/delta + PsiHxY(i,j));
    end
        
    % ========================= By =============================
    % Hy Psi array.
    i=1:IHy-1;
    j=1:JHy;
    PsiHyX(i,j) = (Cmx/delta)*(Ez(i+1,j,n0)-Ez(i,j,n0)) + bmx*PsiHyX(i,j);
    PsiHyX(IHy,j) = (Cmx/delta)*(Ez(1,j,n0)-Ez(IHy,j,n0)) + bmx*PsiHyX(IHy,j);
    % By in normal space.
    j=(1+PMLw):JHy-PMLw;
    By(i,j,nf) = By(i,j,n0) + (Ez(i+1,j,n0) - Ez(i,j,n0)) * dt/delta;
    By(IHy,j,nf) = By(IHy,j,n0) + (Ez(1,j,n0) - Ez(IHy,j,n0)) * dt/delta; % PBC
    if PMLw > 0
        % By in lower PML layer.
        i=1:IHy-1;
        j=1:PMLw;
        By(i,j,nf) = By(i,j,n0) + dt*((1/kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1/delta + PsiHyX(i,j));
        By(IHy,j,nf) = By(IHy,j,n0) + dt*((1/kappmx)*(Ez(1,j,n0) - Ez(IHy,j,n0)) * 1/delta + PsiHyX(IHy,j)); % PBC
        % By in upper PML layer.
        i=1:IHy-1;
        j=JHy-PMLw+1:JHy;
        By(i,j,nf) = By(i,j,n0) + dt*((1/kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1/delta + PsiHyX(i,j));
        By(IHy,j,nf) = By(IHy,j,n0) + dt*((1/kappmx)*(Ez(1,j,n0) - Ez(IHy,j,n0)) * 1/delta + PsiHyX(IHy,j)); % PBC
    end    
    
    % --- Sync ---
    
    % BxAve.
    i=3:IHy-3;
    j=(3+PMLw):JHy-PMLw-2;
    BxAve(i,j,nf) = (Bx(i,j,nf)+Bx(i+1,j,nf)+Bx(i,j+1,nf)+Bx(i+1,j+1,nf))/4;
    % ByAve.    
    i=3:IHx-2;
    j=(4+PMLw):(JHx-PMLw-3);
    ByAve(i,j,nf) = (By(i,j,nf)+By(i-1,j,nf)+By(i,j-1,nf)+By(i-1,j-1,nf))/4;
    
    % ========================= Hx =============================
    i=1:IHx;
    % Hx in normal space.
    j=(2+PMLw):(JHx-PMLw-1);
    %Hx(i,j,nf) = Bx(i,j,nf)./u0;
    Hx(i,j,nf) = (ax(i,j,1).*Bx(i,j,nf)+ax(i,j,2).*Bx(i,j,n0)+ax(i,j,3).*Bx(i,j,np)+ax(i,j,4).*ByAve(i,j,nf)+ax(i,j,5).*ByAve(i,j,n0)+ax(i,j,6).*ByAve(i,j,np)-ax(i,j,7).*Hx(i,j,n0)-ax(i,j,8).*Hx(i,j,np))./ax(i,j,9);
    if PMLw > 0
        % Hx in lower PML layer.
        j=2:PMLw+1;
        Hx(i,j,nf) = Bx(i,j,nf)./u0;
        % Hx in upper PML layer.
        j=JHx-PMLw:JHx-1;
        Hx(i,j,nf) = Bx(i,j,nf)./u0;
    end
    
    % ========================= Hy =============================
    % Hy in normal space.
    i=1:IHy;
    j=(1+PMLw):JHy-PMLw;
    %Hy(i,j,nf) = By(i,j,nf)./u0;
    Hy(i,j,nf) = (ay(i,j,1).*By(i,j,nf)+ay(i,j,2).*By(i,j,n0)+ay(i,j,3).*By(i,j,np)+ay(i,j,4).*BxAve(i,j,nf)+ay(i,j,5).*BxAve(i,j,n0)+ay(i,j,6).*BxAve(i,j,np)-ay(i,j,7).*Hy(i,j,n0)-ay(i,j,8).*Hy(i,j,np))./ay(i,j,9);
    if PMLw > 0
        % By in lower PML layer.
        i=1:IHy;
        j=1:PMLw;
        Hy(i,j,nf) = By(i,j,nf)./u0;
        % By in upper PML layer.
        i=1:IHy;
        j=JHy-PMLw+1:JHy;
        Hy(i,j,nf) = By(i,j,nf)./u0;
    end
    
    % ========================= Dz and Ez =============================
    % Psi arrays.
    i=2:IEz;
    j=1:JEz;
    PsiEzX(i,j) = (Cex/delta)*(Hy(i,j,nf)-Hy(i-1,j,nf)) + bex*PsiEzX(i,j);
    PsiEzX(1,j) = (Cex/delta)*(Hy(1,j,nf)-Hy(IEz,j,nf)) + bex*PsiEzX(1,j); % PBC
    PsiEzY(i,j) = (Cey/delta)*(-Hx(i,j+1,nf)+Hx(i,j,nf)) + bey*PsiEzY(i,j);
    PsiEzY(1,j) = (Cey/delta)*(-Hx(1,j+1,nf)+Hx(1,j,nf)) + bey*PsiEzY(1,j); % PBC
    % Dz in Normal Space.
    j=(1+PMLw):(JEz-PMLw);
    Dz(i,j,nf) = Dz(i,j,n0) + (Hy(i,j,nf)-Hy(i-1,j,nf)-Hx(i,j+1,nf)+Hx(i,j,nf)) * dt/delta;
    Dz(1,j,nf) = Dz(1,j,n0) + (Hy(1,j,nf)-Hy(IEz,j,nf)-Hx(1,j+1,nf)+Hx(1,j,nf)) * dt/delta; % PBC
    i=1:IEz;
    Ez(i,j,nf) = az(i,j,1).*(Dz(i,j,nf)-2*Dz(i,j,n0)+Dz(i,j,np))+az(i,j,2).*(Dz(i,j,nf)-Dz(i,j,np))+az(i,j,3).*(2*Ez(i,j,n0)-Ez(i,j,np))+az(i,j,4).*(2*Ez(i,j,n0)+Ez(i,j,np))+az(i,j,5).*Ez(i,j,np);
    if PMLw > 0
        % Dz in lower PML layer.
        i=2:IEz;
        j=1:PMLw;
        Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1/kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1/kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1/delta + PsiEzX(i,j) + PsiEzY(i,j));
        Dz(1,j,nf) = Dz(1,j,n0) + dt*(((1/kappex)*(Hy(1,j,nf)-Hy(IEz,j,nf))+(1/kappey)*(-Hx(1,j+1,nf)+Hx(1,j,nf))) * 1/delta + PsiEzX(1,j) + PsiEzY(1,j)); % PBC
        i=1:IEz;
        Ez(i,j,nf) = Dz(i,j,nf)./e0;
        % Dz in upper PML layer.
        i=2:IEz;
        j=JEz-PMLw+1:JEz;
        Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1/kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1/kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1/delta + PsiEzX(i,j) + PsiEzY(i,j));
        Dz(1,j,nf) = Dz(1,j,n0) + dt*(((1/kappex)*(Hy(1,j,nf)-Hy(IEz,j,nf))+(1/kappey)*(-Hx(1,j+1,nf)+Hx(1,j,nf))) * 1/delta + PsiEzX(1,j) + PsiEzY(1,j)); % PBC
        i=1:IEz;
        Ez(i,j,nf) = Dz(i,j,nf)./e0;
    end
            
    % ====================== Source ===================
    if SourcePlane == 1
        i = 1:IEz;
        j = SourceLocationY;        
    else
        i = SourceLocationX;
        j = SourceLocationY;
    end
    % Source.
    if SourceChoice == 1
        Ez(i,j,nf) = Ez(i,j,nf) + exp( -1*((n-td)/(PulseWidth/4))^2 ) * Sc;
    elseif SourceChoice == 2
        Ez(i,j,nf) = Ez(i,j,nf) + sin(2*pi*f*(n)*dt) * Sc;
    elseif SourceChoice == 3
        Ez(i,j,nf) = Ez(i,j,nf) + (1-2*(pi*fp*(n*dt-dr))^2)*exp(-1*(pi*fp*(n*dt-dr))^2) * Sc;
    end
    Dz(i,j,nf) = e0*Ez(i,j,nf);
    Ez(:,:,nf) = Ez(:,:,nf).*EzMask;
    Dz(:,:,nf) = Dz(:,:,nf).*EzMask;
    % Transmitted fields.
    Ezt(n+1) = Ez(IEz/2,round((2*PMLw+J)/3),nf);
    Eztt(n+1) = Ez(IEz/2,round((2*PMLw+J)/3)+10,nf);
    
    % Fields for calculation of refractive index.
    Ezy1(n+1) = Ez(IEz/2,Y1,nf);
    Ezy2(n+1) = Ez(IEz/2,Y2, nf);
    
    if (SaveFields == 1 && mod(n, SnapshotInterval) == 0)
        EzSnapshots(:,:,n/SnapshotInterval+1) = Ez(1+(0:SnapshotResolution:(IEz-1)), 1+(0:SnapshotResolution:(JEz-1)), nf);
    end
    
    np = mod(np, 3)+1;
    n0 = mod(n0, 3)+1;
    nf = mod(nf, 3)+1;
end
\end{lstlisting}
\section{C++ Simulation Code}
\label{App:C++-Simulation-Cloak}
\lstset{language=[ISO]C++, commentstyle=\color{green!50!black}, keywordstyle=\color{blue}, stringstyle=\color{red!60!black}}
\begin{lstlisting}
int CFDTD2DCloak::RunSimulationCPU(bool SaveFields)
{
	stringstream framestream;
	string basename = "FieldData/Ez";
	string filename;
	fstream snapshot;
	frame = 0U;

	cout << "Simulation (CPU) started..." << endl;
	for (unsigned int n=0; n<MaxTime; n++)
	{
		if (n%SnapshotInterval == 0)
			cout << "\r" << setprecision(4) << (float)n*100/(MaxTime-1) << "%  " << flush;

		// ========================== Bx ==========================
		for (unsigned int i=0; i<IHx; i++)
		{
			// Calculation of PsiHxY.
			for (unsigned int j=1; j<JHx-1; j++)
				PsiHxY(i,j) = (Cmy/delta)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) + bmy*PsiHxY(i,j);

			// Bx in normal space.
			for (unsigned int j=1+PMLw; j<JHx-PMLw-1; j++)
				Bx(i,j,nf) = Bx(i,j,n0) + (-Ez(i,j,n0) + Ez(i,j-1,n0)) * dt/delta;

			// Bx in lower PML.
			for (unsigned int j=1; j<PMLw+1; j++)
				Bx(i,j,nf) = Bx(i,j,n0) + dt*((1./kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1./delta + PsiHxY(i,j));

			// Bx in upper PML.
			for (unsigned int j=JHx-PMLw-1; j<JHx-1; j++)
				Bx(i,j,nf) = Bx(i,j,n0) + dt*((1./kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1./delta + PsiHxY(i,j));
		}

		// ========================== By ==========================
		for (unsigned int i=0; i<IHy-1; i++)
		{
			// PsiHyX arrays.
			for (unsigned int j=0; j<JHy; j++)
			{
				PsiHyX(i,j) = (Cmx/delta)*(Ez(i+1,j,n0)-Ez(i,j,n0)) + bmx*PsiHyX(i,j);
				if (i==0)
					PsiHyX(IHy-1,j) = (Cmx/delta)*(Ez(0,j,n0)-Ez(IHy-1,j,n0)) + bmx*PsiHyX(IHy-1,j); // PBC
			}
			// By in normal space.
			for (unsigned int j=PMLw; j<JHy-PMLw; j++)
			{
				By(i,j,nf) = By(i,j,n0) + (Ez(i+1,j,n0) - Ez(i,j,n0)) * dt/delta;
				if (i==0)
					By(IHy-1,j,nf) = By(IHy-1,j,n0) + (Ez(0,j,n0) - Ez(IHy-1,j,n0)) * dt/delta; // PBC
			}
			// By in Lower PML.
			for (unsigned int j=0; j<PMLw; j++)
			{
				By(i,j,nf) = By(i,j,n0) + dt*((1./kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1./delta + PsiHyX(i,j));
				if (i==0)
					By(IHy-1,j,nf) = By(IHy-1,j,n0) + dt*((1./kappmx)*(Ez(0,j,n0) - Ez(IHy-1,j,n0)) * 1./delta + PsiHyX(IHy-1,j)); // PBC
			}
			// By in upper PML.
			for (unsigned int j=JHy-PMLw; j<JHy; j++)
			{
				By(i,j,nf) = By(i,j,n0) + dt*((1./kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1./delta + PsiHyX(i,j));
				if (i==0)
					By(IHy-1,j,nf) = By(IHy-1,j,n0) + dt*((1./kappmx)*(Ez(0,j,n0) - Ez(IHy-1,j,n0)) * 1./delta + PsiHyX(IHy-1,j)); // PBC
			}
		}

		// --- Sync ---

		for (unsigned int i=2; i<IHy-2; i++)
		{
			for (unsigned int j=3+PMLw; j<JHx-PMLw-3; j++)
			{
				BxAve(i,j,nf) = (Bx(i,j,nf)+Bx(i+1,j,nf)+Bx(i,j+1,nf)+Bx(i+1,j+1,nf))/4.;
				ByAve(i,j,nf) = (By(i,j,nf)+By(i-1,j,nf)+By(i,j-1,nf)+By(i-1,j-1,nf))/4.;
			}
		}

		// ========================== Hx ==========================
		for (unsigned int i=0; i<IHx; i++)
		{
			// Hx in normal space.
			for (unsigned int j=1+PMLw; j<JHx-PMLw-1; j++)
			{
				//Hx(i,j,nf) = Bx(i,j,nf)/u0;
				Hx(i,j,nf) = (ax(i,j,0)*Bx(i,j,nf)+ax(i,j,1)*Bx(i,j,n0)+ax(i,j,2)*Bx(i,j,np)+ax(i,j,3)*ByAve(i,j,nf)+ax(i,j,4)*ByAve(i,j,n0)+ax(i,j,5)*ByAve(i,j,np)-ax(i,j,6)*Hx(i,j,n0)-ax(i,j,7)*Hx(i,j,np))/ax(i,j,8);
			}
			// Hx in lower PML.
			for (unsigned int j=1; j<PMLw+1; j++)
				Hx(i,j,nf) = Bx(i,j,nf)/u0;

			// Hx in upper PML.
			for (unsigned int j=JHx-PMLw-1; j<JHx-1; j++)
				Hx(i,j,nf) = Bx(i,j,nf)/u0;
		}

		// ========================== Hy ==========================
		for (unsigned int i=0; i<IHy-1; i++)
		{
			// Hy in normal space.
			for (unsigned int j=PMLw; j<JHy-PMLw; j++)
			{
				//Hy(i,j,nf) = By(i,j,nf)/u0;
				Hy(i,j,nf) = (ay(i,j,0)*By(i,j,nf)+ay(i,j,1)*By(i,j,n0)+ay(i,j,2)*By(i,j,np)+ay(i,j,3)*BxAve(i,j,nf)+ay(i,j,4)*BxAve(i,j,n0)+ay(i,j,5)*BxAve(i,j,np)-ay(i,j,6)*Hy(i,j,n0)-ay(i,j,7)*Hy(i,j,np))/ay(i,j,8);
			}
			// Hy in Lower PML.
			for (unsigned int j=0; j<PMLw; j++)
			{
				Hy(i,j,nf) = By(i,j,nf)/u0;
				if (i==0)
					Hy(IHy-1,j,nf) = By(IHy-1,j,nf)/u0; // PBC
			}
			// Hy in upper PML.
			for (unsigned int j=JHy-PMLw; j<JHy; j++)
			{
				Hy(i,j,nf) = By(i,j,nf)/u0;
				if (i==0)
					Hy(IHy-1,j,nf) = By(IHy-1,j,nf)/u0; // PBC
			}
		}

		// ========================== Dz and Ez ==========================
		for (unsigned int i=1; i<IEz; i++)
		{
			// Psi arrays.
			for (unsigned int j=0; j<JEz; j++)
			{
				PsiEzX(i,j) = (Cex/delta)*(Hy(i,j,nf)-Hy(i-1,j,nf)) + bex*PsiEzX(i,j);
				PsiEzY(i,j) = (Cey/delta)*(-Hx(i,j+1,nf)+Hx(i,j,nf)) + bey*PsiEzY(i,j);
				if (i==1)
				{
					PsiEzX(0,j) = (Cex/delta)*(Hy(0,j,nf)-Hy(IEz-1,j,nf)) + bex*PsiEzX(0,j); // PBC
					PsiEzY(0,j) = (Cey/delta)*(-Hx(0,j+1,nf)+Hx(0,j,nf)) + bey*PsiEzY(0,j); // PBC
				}
			}
			// Dz in normal space.
			for (unsigned int j=PMLw; j<JEz-PMLw; j++)
			{
				Dz(i,j,nf) = Dz(i,j,n0) + (Hy(i,j,nf)-Hy(i-1,j,nf)-Hx(i,j+1,nf)+Hx(i,j,nf)) * dt/delta;
				Ez(i,j,nf) = EzMask(i,j) * (az(i,j,0)*(Dz(i,j,nf)-2.*Dz(i,j,n0)+Dz(i,j,np))+az(i,j,1)*(Dz(i,j,nf)-Dz(i,j,np))+az(i,j,2)*(2.*Ez(i,j,n0)-Ez(i,j,np))+az(i,j,3)*(2.*Ez(i,j,n0)+Ez(i,j,np))+az(i,j,4)*Ez(i,j,np));
				if (i==1)
				{
					Dz(0,j,nf) = Dz(0,j,n0) + (Hy(0,j,nf)-Hy(IEz-1,j,nf)-Hx(0,j+1,nf)+Hx(0,j,nf)) * dt/delta; // PBC
					Ez(0,j,nf) = EzMask(0,j) * (az(0,j,0)*(Dz(0,j,nf)-2.*Dz(0,j,n0)+Dz(0,j,np))+az(0,j,1)*(Dz(0,j,nf)-Dz(0,j,np))+az(0,j,2)*(2.*Ez(0,j,n0)-Ez(0,j,np))+az(0,j,4)*(2.*Ez(0,j,n0)+Ez(0,j,np))+az(0,j,4)*Ez(0,j,np)); // PBC
				}
			}
			// Dz in lower PML.
			for (unsigned int j=0; j<PMLw; j++)
			{
				Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1./kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1./kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1./delta + PsiEzX(i,j) + PsiEzY(i,j));
				Ez(i,j,nf) = EzMask(i,j) * Dz(i,j,nf)/e0;
				if (i==1)
				{
					Dz(0,j,nf) = Dz(0,j,n0) + dt*(((1./kappex)*(Hy(0,j,nf)-Hy(IEz-1,j,nf))+(1./kappey)*(-Hx(0,j+1,nf)+Hx(0,j,nf))) * 1./delta + PsiEzX(0,j) + PsiEzY(0,j)); // PBC
					Ez(0,j,nf) = EzMask(0,j) * Dz(0,j,nf)/e0; // PBC
				}
			}
			// Dz in upper PML.
			for (unsigned int j=JEz-PMLw; j<JEz; j++)
			{
				Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1./kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1./kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1./delta + PsiEzX(i,j) + PsiEzY(i,j));
				Ez(i,j,nf) = EzMask(0,j) * Dz(i,j,nf)/e0;
				if (i==1)
				{
					Dz(0,j,nf) = Dz(0,j,n0) + dt*(((1./kappex)*(Hy(0,j,nf)-Hy(IEz-1,j,nf))+(1./kappey)*(-Hx(0,j+1,nf)+Hx(0,j,nf))) * 1./delta + PsiEzX(0,j) + PsiEzY(0,j)); // PBC
					Ez(0,j,nf) = EzMask(0,j) * Dz(0,j,nf)/e0; // PBC
				}
			}
		}

		// ====================== Source ======================
		if (SourcePlane == 1)
		{
			unsigned int j=SourceLocationY;
			for (unsigned int i=0; i<IEz; i++)
			{
				if (SourceChoice == 1)
					Ez(i,j,nf) += exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
				else if (SourceChoice == 2)
					Ez(i,j,nf) += sin(2.*pi*f*(PRECISION)n*dt) * Sc;
				else if (SourceChoice == 3)
					Ez(i,j,nf) += (1.-2.*pow(pi*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(pi*fp*((PRECISION)n*dt-dr),2)) * Sc;

				Dz(i,j,nf) = e0*Ez(i,j,nf);
			}
		}
		else
		{
			unsigned int i=SourceLocationX;
			unsigned int j=SourceLocationY;
			if (SourceChoice == 1)
				Ez(i,j,nf) += exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
			else if (SourceChoice == 2)
				Ez(i,j,nf) += sin(2.*pi*f*(PRECISION)n*dt) * Sc;
			else if (SourceChoice == 3)
				Ez(i,j,nf) += (1.-2.*pow(pi*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(pi*fp*((PRECISION)n*dt-dr),2)) * Sc;

			Dz(i,j,nf) = e0*Ez(i,j,nf);
		}
		// Recording transmitted fields.
		Ezt[n] = Ez(IEz/2,x1,nf);
		Eztt[n] = Ez(IEz/2,x1+1,nf);
		// Fields for refractive index.
		Ezy1[n] = Ez(IEz/2,Y1,nf);
		Ezy2[n] = Ez(IEz/2,Y2,nf);

		// Saving electric field snapshot.
		if (n%SnapshotInterval == 0 && SaveFields == true)
		{
			// Write E-field to file.
			framestream.str(std::string());			// Clearing stringstream contents.
			framestream << ++frame;
			filename = basename + framestream.str() + ".fdt";
			snapshot.open(filename.c_str(), std::ios::out|std::ios::binary);
			snapshot.write((char*)&(Ez(0,0,nf)), sizeof(PRECISION)*IEz*JEz);
			snapshot.close();
		}

		np = (np+1)%3;
		n0 = (n0+1)%3;
		nf = (nf+1)%3;
	}
	// Saving electric field snapshot.
	if (SaveFields == true)
	{
		fstream parametersfile;
		parametersfile.open("FieldData/Parameters.smp", std::ios::out|std::ios::binary|std::ios::app);
		parametersfile.write((char*)&(frame), sizeof(unsigned int));
		parametersfile.close();
		// Write saved fields to files.
		snapshot.open("FieldData/Ezi.fdt", std::ios::out|std::ios::binary);
		snapshot.write((char*)Ezi, sizeof(PRECISION)*MaxTime);
		snapshot.close();
		snapshot.open("FieldData/Ezt.fdt", std::ios::out|std::ios::binary);
		snapshot.write((char*)Ezt, sizeof(PRECISION)*MaxTime);
		snapshot.close();
		snapshot.open("FieldData/Eztt.fdt", std::ios::out|std::ios::binary);
		snapshot.write((char*)Eztt, sizeof(PRECISION)*MaxTime);
		snapshot.close();
		snapshot.open("FieldData/Ezy1.fdt", std::ios::out|std::ios::binary);
		snapshot.write((char*)Ezy1, sizeof(PRECISION)*MaxTime);
		snapshot.close();
		snapshot.open("FieldData/Ezy2.fdt", std::ios::out|std::ios::binary);
		snapshot.write((char*)Ezy2, sizeof(PRECISION)*MaxTime);
		snapshot.close();
	}
	cout << endl << "Simulation (CPU) completed!" << endl;
	return 0;
}
\end{lstlisting}
\section{CUDA Simulation Kernels Code}
\label{App:CUDA-Cloak-Simulation-Kernels}
\begin{lstlisting}
// Simulations kernels.
template <unsigned int BlockX, unsigned int BlockY> __global__ void FDTD2DCloakKernel_Simulation_B(
				const unsigned int I,
				const unsigned int J,
				const unsigned int PMLw,
				const unsigned int PulseWidth,
				const unsigned int td,
				const unsigned int SourceChoice,
				const unsigned int SourcePlane,
				const unsigned int SourceLocationX,
				const unsigned int SourceLocationY,
				const PRECISION c,
				const PRECISION pi,
				const PRECISION e0,
				const PRECISION u0,
				const PRECISION dt,
				const PRECISION delta,
				const PRECISION Sc,
				const PRECISION f,
				const PRECISION fp,
				const PRECISION dr,
				const unsigned int IEz, const unsigned int JEz,
				const unsigned int IHx, const unsigned int JHx,
				const unsigned int IHy, const unsigned int JHy,
				const PRECISION *Ez_, const PRECISION *Dz_, PRECISION *Hx_, PRECISION *Bx_, PRECISION *Hy_, PRECISION *By_,
				const PRECISION *ax_, const PRECISION *ay_, const PRECISION *az_,
				const PRECISION *PsiEzX_, const PRECISION *PsiEzY_, PRECISION *PsiHyX_, PRECISION *PsiHxY_,
				const PRECISION kappex, const PRECISION kappey, const PRECISION kappmx, const PRECISION kappmy,
				const PRECISION bex, const PRECISION bey, const PRECISION bmx, const PRECISION bmy,
				const PRECISION Cex, const PRECISION Cey, const PRECISION Cmx, const PRECISION Cmy,
				PRECISION *Ezt, PRECISION *Eztt, PRECISION *Ezy1, PRECISION *Ezy2,
				const unsigned int x1, const unsigned int Y1, const unsigned int Y2,
				const unsigned int n,
				const unsigned int np,
				const unsigned int n0,
				const unsigned int nf)
{
	const unsigned int i = BlockX*blockIdx.x+threadIdx.x;
	const unsigned int j = BlockY*blockIdx.y+threadIdx.y;

	// ========================== Bx ==========================
	// Calculation of PsiHxY.
	if (j>0 && j<JHx-1)
		PsiHxY(i,j) = (Cmy/delta)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) + bmy*PsiHxY(i,j);

	// Bx in normal space.
	if (j>PMLw && j<JHx-PMLw-1)
		Bx(i,j,nf) = Bx(i,j,n0) + (-Ez(i,j,n0) + Ez(i,j-1,n0)) * dt/delta;

		// Bx in lower PML.
	if (j>0 && j<PMLw+1)
		Bx(i,j,nf) = Bx(i,j,n0) + dt*((1./kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1./delta + PsiHxY(i,j));

	// Bx in upper PML.
	if (j>JHx-PMLw-2 && j<JHx-1)
		Bx(i,j,nf) = Bx(i,j,n0) + dt*((1./kappmy)*(-Ez(i,j,n0) + Ez(i,j-1,n0)) * 1./delta + PsiHxY(i,j));

	// ========================== By ==========================
	if (i<IHy-1)
	{
		// PsiHyX array.
		if (j<JHy)
			PsiHyX(i,j) = (Cmx/delta)*(Ez(i+1,j,n0)-Ez(i,j,n0)) + bmx*PsiHyX(i,j);

		// By in normal space.
		if (j>PMLw-1 && j<JHy-PMLw)
			By(i,j,nf) = By(i,j,n0) + (Ez(i+1,j,n0) - Ez(i,j,n0)) * dt/delta;

			// By in Lower PML.
		if (j<PMLw)
			By(i,j,nf) = By(i,j,n0) + dt*((1./kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1./delta + PsiHyX(i,j));

		// By in upper PML.
		if (j>JHy-PMLw-1 && j<JHy)
			By(i,j,nf) = By(i,j,n0) + dt*((1./kappmx)*(Ez(i+1,j,n0) - Ez(i,j,n0)) * 1./delta + PsiHyX(i,j));
	}
	else
	{
		// PsiHyX array.
		if (j<JHy)
			PsiHyX(IHy-1,j) = (Cmx/delta)*(Ez(0,j,n0)-Ez(IHy-1,j,n0)) + bmx*PsiHyX(IHy-1,j); // PBC

		// By in normal space.
		if (j>PMLw-1 && j<JHy-PMLw)
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + (Ez(0,j,n0) - Ez(IHy-1,j,n0)) * dt/delta; // PBC

		// By in Lower PML.
		if (j<PMLw)
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + dt*((1./kappmx)*(Ez(0,j,n0) - Ez(IHy-1,j,n0)) * 1./delta + PsiHyX(IHy-1,j)); // PBC

		// By in upper PML.
		if (j>JHy-PMLw-1 && j<JHy)
			By(IHy-1,j,nf) = By(IHy-1,j,n0) + dt*((1./kappmx)*(Ez(0,j,n0) - Ez(IHy-1,j,n0)) * 1./delta + PsiHyX(IHy-1,j)); // PBC
	}
}
template <unsigned int BlockX, unsigned int BlockY> __global__ void FDTD2DCloakKernel_Simulation_H(
				const unsigned int I,
				const unsigned int J,
				const unsigned int PMLw,
				const unsigned int PulseWidth,
				const unsigned int td,
				const unsigned int SourceChoice,
				const unsigned int SourcePlane,
				const unsigned int SourceLocationX,
				const unsigned int SourceLocationY,
				const PRECISION c,
				const PRECISION pi,
				const PRECISION e0,
				const PRECISION u0,
				const PRECISION dt,
				const PRECISION delta,
				const PRECISION Sc,
				const PRECISION f,
				const PRECISION fp,
				const PRECISION dr,
				const unsigned int IEz, const unsigned int JEz,
				const unsigned int IHx, const unsigned int JHx,
				const unsigned int IHy, const unsigned int JHy,
				const PRECISION *Ez_, const PRECISION *Dz_, PRECISION *Hx_, const PRECISION *Bx_, PRECISION *BxAve_, PRECISION *Hy_, const PRECISION *By_, PRECISION *ByAve_,
				const PRECISION *ax_, const PRECISION *ay_, const PRECISION *az_,
				const PRECISION *PsiEzX_, const PRECISION *PsiEzY_, PRECISION *PsiHyX_, PRECISION *PsiHxY_,
				const PRECISION kappex, const PRECISION kappey, const PRECISION kappmx, const PRECISION kappmy,
				const PRECISION bex, const PRECISION bey, const PRECISION bmx, const PRECISION bmy,
				const PRECISION Cex, const PRECISION Cey, const PRECISION Cmx, const PRECISION Cmy,
				PRECISION *Ezt, PRECISION *Eztt, PRECISION *Ezy1, PRECISION *Ezy2,
				const unsigned int x1, const unsigned int Y1, const unsigned int Y2,
				const unsigned int n,
				const unsigned int np,
				const unsigned int n0,
				const unsigned int nf)
{
	const unsigned int i = BlockX*blockIdx.x+threadIdx.x;
	const unsigned int j = BlockY*blockIdx.y+threadIdx.y;

	// ========================== BxAve and ByAve ==========================
	if (i>1 && i<IHy-2 && j>2+PMLw && j<JHx-PMLw-3)
	{
		BxAve(i,j,nf) = (Bx(i,j,nf)+Bx(i+1,j,nf)+Bx(i,j+1,nf)+Bx(i+1,j+1,nf))/4.;
		ByAve(i,j,nf) = (By(i,j,nf)+By(i-1,j,nf)+By(i,j-1,nf)+By(i-1,j-1,nf))/4.;
	}

	// ========================== Hx ==========================
	// Hx in normal space.
	if (j>PMLw && j<JHx-PMLw-1)
		Hx(i,j,nf) = (ax(i,j,0)*Bx(i,j,nf)+ax(i,j,1)*Bx(i,j,n0)+ax(i,j,2)*Bx(i,j,np)+ax(i,j,3)*ByAve(i,j,nf)+ax(i,j,4)*ByAve(i,j,n0)+ax(i,j,5)*ByAve(i,j,np)-ax(i,j,6)*Hx(i,j,n0)-ax(i,j,7)*Hx(i,j,np))/ax(i,j,8);

	// Bx in lower PML.
	if (j>0 && j<PMLw+1)
		Hx(i,j,nf) = Bx(i,j,nf)/u0;

	// Bx in upper PML.
	if (j>JHx-PMLw-2 && j<JHx-1)
		Hx(i,j,nf) = Bx(i,j,nf)/u0;

	// ========================== Hy ==========================
	// Hy in normal space.
	if (j>PMLw-1 && j<JHy-PMLw)
		Hy(i,j,nf) = (ay(i,j,0)*By(i,j,nf)+ay(i,j,1)*By(i,j,n0)+ay(i,j,2)*By(i,j,np)+ay(i,j,3)*BxAve(i,j,nf)+ay(i,j,4)*BxAve(i,j,n0)+ay(i,j,5)*BxAve(i,j,np)-ay(i,j,6)*Hy(i,j,n0)-ay(i,j,7)*Hy(i,j,np))/ay(i,j,8);

	// Hy in Lower PML.
	if (j<PMLw)
		Hy(i,j,nf) = By(i,j,nf)/u0;

	// Hy in upper PML.
	if (j>JHy-PMLw-1 && j<JHy)
		Hy(i,j,nf) = By(i,j,nf)/u0;
}
template <unsigned int BlockX, unsigned int BlockY> __global__ void FDTD2DCloakKernel_Simulation_E(
				const unsigned int I,
				const unsigned int J,
				const unsigned int PMLw,
				const unsigned int PulseWidth,
				const unsigned int td,
				const unsigned int SourceChoice,
				const unsigned int SourcePlane,
				const unsigned int SourceLocationX,
				const unsigned int SourceLocationY,
				const PRECISION c,
				const PRECISION pi,
				const PRECISION e0,
				const PRECISION u0,
				const PRECISION dt,
				const PRECISION delta,
				const PRECISION Sc,
				const PRECISION f,
				const PRECISION fp,
				const PRECISION dr,
				const unsigned int IEz, const unsigned int JEz,
				const unsigned int IHx, const unsigned int JHx,
				const unsigned int IHy, const unsigned int JHy,
				PRECISION *Ez_, PRECISION *Dz_, const PRECISION *EzMask_, const PRECISION *Hx_, const PRECISION *Bx_, const PRECISION *Hy_, const PRECISION *By_,
				const PRECISION *ax_, const PRECISION *ay_, const PRECISION *az_,
				PRECISION *PsiEzX_, PRECISION *PsiEzY_, const PRECISION *PsiHyX_, const PRECISION *PsiHxY_,
				const PRECISION kappex, const PRECISION kappey, const PRECISION kappmx, const PRECISION kappmy,
				const PRECISION bex, const PRECISION bey, const PRECISION bmx, const PRECISION bmy,
				const PRECISION Cex, const PRECISION Cey, const PRECISION Cmx, const PRECISION Cmy,
				PRECISION *Ezt, PRECISION *Eztt, PRECISION *Ezy1, PRECISION *Ezy2,
				const unsigned int x1, const unsigned int Y1, const unsigned int Y2,
				const unsigned int n,
				const unsigned int np,
				const unsigned int n0,
				const unsigned int nf)
{
	const unsigned int i = BlockX*blockIdx.x+threadIdx.x;
	const unsigned int j = BlockY*blockIdx.y+threadIdx.y;

	// ========================== Dz and Ez ==========================
	if (i>0)
	{
		// Psi arrays.
		if (j<JEz)
		{
			PsiEzX(i,j) = (Cex/delta)*(Hy(i,j,nf)-Hy(i-1,j,nf)) + bex*PsiEzX(i,j);
			PsiEzY(i,j) = (Cey/delta)*(-Hx(i,j+1,nf)+Hx(i,j,nf)) + bey*PsiEzY(i,j);
		}
		// Dz in normal space.
		if (j>PMLw-1 && j<JEz-PMLw)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + (Hy(i,j,nf)-Hy(i-1,j,nf)-Hx(i,j+1,nf)+Hx(i,j,nf)) * dt/delta;
			Ez(i,j,nf) = EzMask(i,j) * (az(i,j,0)*(Dz(i,j,nf)-2.*Dz(i,j,n0)+Dz(i,j,np))+az(i,j,1)*(Dz(i,j,nf)-Dz(i,j,np))+az(i,j,2)*(2.*Ez(i,j,n0)-Ez(i,j,np))+az(i,j,3)*(2.*Ez(i,j,n0)+Ez(i,j,np))+az(i,j,4)*Ez(i,j,np));
		}
		// Dz in lower PML.
		if (j<PMLw)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1./kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1./kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1./delta + PsiEzX(i,j) + PsiEzY(i,j));
			Ez(i,j,nf) = Dz(i,j,nf)/e0;
		}
		// Dz in upper PML.
		if (j>JEz-PMLw-1 && j<JEz)
		{
			Dz(i,j,nf) = Dz(i,j,n0) + dt*(((1./kappex)*(Hy(i,j,nf)-Hy(i-1,j,nf))+(1./kappey)*(-Hx(i,j+1,nf)+Hx(i,j,nf))) * 1./delta + PsiEzX(i,j) + PsiEzY(i,j));
			Ez(i,j,nf) = Dz(i,j,nf)/e0;
		}
	}
	else
	{
		// Psi arrays.
		if (j<JEz)
		{
			PsiEzX(0,j) = (Cex/delta)*(Hy(0,j,nf)-Hy(IEz-1,j,nf)) + bex*PsiEzX(0,j); // PBC
			PsiEzY(0,j) = (Cey/delta)*(-Hx(0,j+1,nf)+Hx(0,j,nf)) + bey*PsiEzY(0,j); // PBC
		}
		// Dz in normal space.
		if (j>PMLw-1 && j<JEz-PMLw)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + (Hy(0,j,nf)-Hy(IEz-1,j,nf)-Hx(0,j+1,nf)+Hx(0,j,nf)) * dt/delta; // PBC
			Ez(0,j,nf) = EzMask(0,j) * (az(0,j,0)*(Dz(0,j,nf)-2.*Dz(0,j,n0)+Dz(0,j,np))+az(0,j,1)*(Dz(0,j,nf)-Dz(0,j,np))+az(0,j,2)*(2.*Ez(0,j,n0)-Ez(0,j,np))+az(0,j,4)*(2.*Ez(0,j,n0)+Ez(0,j,np))+az(0,j,4)*Ez(0,j,np)); // PBC
		}
		// Dz in lower PML.
		if (j<PMLw)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + dt*(((1./kappex)*(Hy(0,j,nf)-Hy(IEz-1,j,nf))+(1./kappey)*(-Hx(0,j+1,nf)+Hx(0,j,nf))) * 1./delta + PsiEzX(0,j) + PsiEzY(0,j)); // PBC
			Ez(0,j,nf) = Dz(0,j,nf)/e0; // PBC
		}
		// Dz in upper PML.
		if (j>JEz-PMLw-1 && j<JEz)
		{
			Dz(0,j,nf) = Dz(0,j,n0) + dt*(((1./kappex)*(Hy(0,j,nf)-Hy(IEz-1,j,nf))+(1./kappey)*(-Hx(0,j+1,nf)+Hx(0,j,nf))) * 1./delta + PsiEzX(0,j) + PsiEzY(0,j)); // PBC
			Ez(0,j,nf) = Dz(0,j,nf)/e0; // PBC
		}
	}

	// ====================== Source ======================
	if (SourcePlane == 1 && i<IEz && j == SourceLocationY)
	{
		if (SourceChoice == 1)
			Ez(i,j,nf) += exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
		else if (SourceChoice == 2)
			Ez(i,j,nf) += sin(2.*pi*f*(PRECISION)n*dt) * Sc;
		else if (SourceChoice == 3)
			Ez(i,j,nf) += (1.-2.*pow(pi*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(pi*fp*((PRECISION)n*dt-dr),2)) * Sc;

		Dz(i,j,nf) = e0*Ez(i,j,nf);
	}
	else if (i == SourceLocationX && j == SourceLocationY)
	{
		if (SourceChoice == 1)
			Ez(i,j,nf) += exp(-1.*pow(((PRECISION)n-(PRECISION)td)/((PRECISION)PulseWidth/4.),2)) * Sc;
		else if (SourceChoice == 2)
			Ez(i,j,nf) += sin(2.*pi*f*(PRECISION)n*dt) * Sc;
		else if (SourceChoice == 3)
			Ez(i,j,nf) += (1.-2.*pow(pi*fp*((PRECISION)n*dt-dr),2))*exp(-1.*pow(pi*fp*((PRECISION)n*dt-dr),2)) * Sc;

		Dz(i,j,nf) = e0*Ez(i,j,nf);
	}
	// Recording transmitted fields.
	if (i==IEz/2 && j==x1)
		Ezt[n] = Ez(i,j,nf);
	if (i==IEz/2 && j==2*J/3+PMLw+1)
		Eztt[n] = Ez(i,j,nf);
	// Fields for refractive index.
	if (i==IEz/2 && j==Y1)
		Ezy1[n] = Ez(i,j,nf);
	if (i==IEz/2 && j==Y2)
		Ezy2[n] = Ez(i,j,nf);
}
\end{lstlisting}